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Negative Camber Airfoils for Transonic Flight 


THOMAS B. WOERSCHING* 
The Glenn L. Martin Company 


SUMMARY 


The potentialities of negative camber as a means of increasing 
airfoil \/.’s are pointed out and substantiated by the negative 
lift range of available unclassified transonic test data for positive 
camber airfoils. “W” camber for transonic airfoils 
is proposed as a practical refinement of the “S” 
yields the highest possible A/, but suffers from intolerably large 
negative values of C,,. The NACA 847B110 airfoil (when in- 
verted) has the suggested ‘‘W”’ camber line, is the best transonic 
airfoil for which high subsonic test data are available, and is proof 
of the proposed ‘‘W”’ camber concept. 

Negative camber airfoils with positive values of C,,. are pro- 
posed to eliminate longitudinal ‘“‘tuck under” instability at tran- 
Negative camber airfoils provide greater structural 


The use of 
cember which 


sonic /)’s. 
thickness than a thinner symmetrical airfoil having the same upper 
surface curvature, C),, and M,. 

The theoretical characteristics at Wo 
airfoils derived for use at transonic /)’s are presented. 
0.2 


0 of two ‘“‘W” camber 
The esti- 


mated J/.’s are about 0.864 at C; 


SYMBOLS 


airfoil angle of attack 

two-dimensional angle of attack 

airfoil chord 

airfoil section profile drag coefficient 
airfoil section lift coefficient 

airplane lift coefficient 

airfoil design lift coefficient 

maximum airfoil section lift coefficient 
increment of C; due to angle of attack 
airfoil section pitching moment about c/4 
coefficient 


4Cia 
Cm 
Crm 


airfoil section pitching moment 
Ci =0 


airplane pitching moment coefficient at Cr, = 0 


Cu, 
dCy/dCr 
M 

My 

M, 


airplane pitching moment slope 

local flow Mach Number 

free-stream Mach Number 

free-stream critical Mach Number for local shock 
and/or force divergence 

airfoil local static pressure relative to the atmos- 


Ap 
phere 
local pressure coefficient, Ap/g 


P 
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= Pat which local sonic velocity obtains 
dynamic pressure, (1/2) p V? 
airfoil thickness, fraction of chord 


P. 


velocity 
distance from leading edge along airfoil chord 


air mass density 


INTRODUCTION 


SB prego THE RECENT POSTWAR YEARS, the problem 
of airfoil development for transonic flight appears 
to have been de-emphasized in deference to wing sweep- 
back and supersonic research. As many future air- 
craft will undoubtedly be built to cruise at speeds near 
sonic velocity, it is apparent that the development of 
airfoils for transonic flight should be continued. 

One potentiality in airfoil design that has not yet 
been exploited is the use of negative camber and/or 
combined negative and positive camber. The ad- 
vantages of this design concept have been suggested 
by trends of transonic test data for positive camber 
airfoils in the negative lift range. The camber line ad- 
vocated in this paper is one with “W” in 
which the camber of the forward 75 to 80 per cent of 
the chord is negative, the camber at about x/c = 
0.85 is positive, and the camber aft to the trailing edge 


camber, 


O 


is negative. 

The purpose of this paper is to discuss qualitatively 
the physical basis for this design idea and to point out 
its several advantages and, where possible, to verify it 
with the available unclassified test data. Two analyt- 
ically derived airfoil profiles are proposed in the hope 
of interesting agencies having sufficient funds and facili- 
ties in exploring the concept further. In the opinion 
of the writer, merits of the ““W’’ camber line are proved 
by the negative lift range of the test data available for 
an N.A.C.A. airfoil having an ‘“‘M’’-shaped camber 
line.° 

It is not contended here that the need for sweepback 
can be eliminated by ‘“‘W”’ camber airfoils, but, for 
equal aircraft performance, the amount of sweepback 
needed can be reduced with consequent mitigation of 
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the stall and stability problems. Conversely, at con- 
stant sweepback, the aircraft performance can be sub- 
stantially improved by a negative camber airfoil. 


NEGATIVE CAMBER ON AIRFOIL 
CHARACTERISTICS 


THE EFFECTS OF 


Critical Mach Number 


It is well established that the lower or first divergence 
M, of a conventional unsweptback positive camber or 
symmetrical airfoil is caused by the local potential flow 
over the top surface of the airfoil. The induced ve- 
locity over the top surface causes sonic or supersonic 
flow over the airfoil while J) is subsonic. The //, ob- 
tains when the local flow J/ is sufficiently supersonic to 
cause a shock. 

The pressure distributions on the upper and lower 
surfaces of an airfoil can be integrated separately to 
obtain upper and lower surface C;’s. This was done 
with the dive flight data for a 30-in. chord NACA 
65-415 airfoil, presented in Fig. 1.* It will be noted 
that the airfoil /, is determined by the discontinuity 
in the upper surface integrated C, at J) = 0.68, while 
the variation of the integrated lower surface C, is regu- 
lar and continuous. 


JUNE, 1951 

Similarly, a second or higher divergence .V/, is pro- 
duced when the velocities under the lower surface are 
sufficiently supersonic to cause a shock. This differ- 
ence between the upper and lower surface ./, causes a 
double discontinuity to form the familiar “‘S’’ variation 
in unsweptback airfoil transonic lift and pitching mo- 
ment characteristics. 

Inverting the airfoil camber reduces the mid-chord 
local velocities over the airfoil top surface and increases 
those under the bottom surface. This effect tends to 
make the two /,’s approach each other,’ as illustrated 
schematically in Figs. 2 and 3. The effect of the cam- 
ber inversion is to reduce the airfoil lift at a constant 
chord line a, but lift equivalent to that of a symmetrical 
or positive camber airfoil is obtainable at a slightly 
higher a. 

At the same high speed C,, the upper surface .1/. will 
be higher for a negative camber airfoil than for a sym- 
metrical or positive camber airfoil having the same 
thickness distribution because of less upper surface 
This is true in spite of 


) 


curvature, as shown in Fig. 8. 
the fact that at low subsonic J) a high negative P 
occurs at the negative camber airfoil leading edge. 
This causes a premature shock of small intensity over 
a small chordwise extent of the airfoil leading edge 
which has negligible effect on the airfoil lift, drag, and 
pitching moment characteristics. 

These effects are illustrated by the dive flight data’ 
presented in Fig. 4 for a constant 30-in. chord NACA 
65-411 airfoil. Chordwise pressure distribution plots 
at constant a, = —3.0° for My = 0.418, 0.668, 0.747, 
and 0.800 are presented in conjunction with the cor- 
responding variation of C, with /, at constant a. At 


M, = 0.418, the lower surface leading edge P = —2.2, 
while the P, for sonic local velocity® is —3.38; 1.e., the 
local flow is subsonic. However, at j1/) = 0.668, the 


leading-edge lower surface P = —1.90, while P, = 
— 0.92, indicating that the local flow is supersonic over 
the leading edge to x/c = 0.09. Yet, at 1J = 0.668, 
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NEGATIVE CAMBER AIRFOILS FOR TRANSONIC 


the small leading-edge shock has no effect on the air- 
foil characteristics. At the drag divergence WV, = 
0.747, the leading-edge maximum P = —1.55, with 
P, = —0.6. Thus, the supersonic shock flow extends 
aft to about x/c = 0.32 before the drag divergence M/, 
obtains. At the supercritical My = 0.800, P, = 
—(.44, while the test P exceeds this value over the for- 
ward 49 per cent of the airfoil chord. These data illus- 
trate that a considerable extent of leading-edge shock 
can be tolerated with no adverse effect on the airfoil 
characteristics. 

Rather, it is the shock that usually occurs near the 
airfoil mid-chord at W/ greater than 0.6 for C;’s up to 
about 0.8 which causes the familiar transonic discon- 
tinuities in C,; and C,, and the divergence in Cy as a 
function of 17. However, as C, increases, the leading- 
edge upper surface shock intensity and extent increase 
until, finally, it affects the entire upper surface flow 
and causes the familiar discontinuity in the variation 
of divergence M, with C), indicated in Fig. 8. It is 
logical to expect that an airfoil could be derived which 
would have the maximum negative P for the upper and 
lower surface equal in magnitude and at the same 
chordwise station for one design C;. Such an airfoil 
would have equal upper surface and lower surface 
M,’s, with the shocks occurring at the same chordwise 
station such that the upper and lower surface pressure 
variations should tend to compensate to yield essen- 
tially constant values of C, and C,, for the complete 
airfoil through the transonic range of Mp. 

The effect of negative camber on the variation of 
C, with My should be to increase the C, divergence 1/, 
by the amount of the upper surface increment in J/,. 
However, the slope of the supercritical drag rise should 
be equal to that of a conventional positive camber 
airfoil, and the drag divergence should be more abrupt 
because upper surface and lower surface shocks would 
occur at essentially the same J/». 


Zero Lift Pitching Moment 


A negative camber airfoil would contribute an incre- 
ment of the positive Cy,, required to trim every stable 
airplane configuration at a positive C,. In a conven- 
tional airplane, the + Cy, is provided by the stabilizer, 
but in a tailless flying wing it must be provided by the 
negative camber resulting from upward deflection of 
wing trailing-edge flaps. A negative camber airfoil 
should provide trimmed pitch equilibrium at less cost 
in drag than deflected elevons, because the positive 
Cy, is obtained from negative camber over a greater 
extent of airfoil chord than that of a deflected trailing- 
edge flap. A properly designed negative-camber tail- 
less airplane configuration would be trimmed at cruise 
C,’s with zero elevon deflection. 


Transonic Longitudinal Stability and Control 


One of the more critical problems of transonic flight 
is the tendency toward catastrophic instability, familiar 


FLIGHT 363 








_ MOMENT TO BE TRIMMED 








Cu S 6 Cn 
+ 7 + 
8 
9 











M=.95 
Positive camber airfoil pitching moments. Fic. 6 


Fic. 5 (left). 
Negative camber airfoil pitching moments. 


(right). 


as “tuck under,” in which a nose-down pitching mo- 
ment is produced as the flight M/») approaches 1.0. 
This phenomenon is produced by the facts that the 
pitching moment dCy,/dC,, becomes more 
negative with increasing subsonic My) and that the 
Cy,, for airplanes with positive camber airfoils, also 
becomes more negative with increasing Mp. 

The combination of these two effects is illustrated 
schematically for the positive camber airfoil by the 
family of typical curves of Cy vs. C, at constant Mo 
in Fig. 5. A representative variation of the moment 
to be trimmed, as a function of C, and Mp, is super- 
imposed. It will be seen that the variations of Cy,, and 
dCy,/dC;, as Mo increases cause the trimmed pitching 
moment curve slope to become more positive, making 
the airplane statically unstable. 

It is logical to expect that, if the camber of the air- 
foil were inverted, the variation of Cy, and dCy,/dC, 
with MM, would be mutually compensating. This effect 
of negative camber is illustrated by the pitching moment 
curves shown in Fig. 6. In the negative C, region, the 
slope, dCy/dC,, of a positive camber airfoil increases 
negatively with increasing subsonic Wo, just as it does 
for positive C,’s, as in Fig. 11. However, the positive 
Cy, increases positively as Mo increases, so that, as the 
pitching moment curves of Cy vs. C, at constant Mo 
are rotated negatively with increasing Jy, the Cy, 
intercepts increase positively. Thus, the two vari- 
ations of dC),/dC, and Cy, with /, tend to compensate 
mutually to produce an essentially linear and stable 
variation of trimmed longitudinal pitching moment as 
a function of C, and Mo. Hence, it is suggested by the 
characteristic trends of transonic airfoil pitching mo- 
ment data that longitudinal ‘‘tuck under’ instability 
can be eliminated by the proper application of negative 
camber in the basic airfoil design. 


slope, 


Maximum Lift at Low Subsonic Mach Numbers 


Generally, negative camber in an airfoil decreases 
its C,,,at low subsonic Mo. This effect results from the 
tendency for negative camber to produce a sharper top 
surface leading edge and a relatively flat upper surface, 
which induces separation because of the relatively ad- 
verse pressure gradient compared to a positive camber 
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Fic. 7. Camber lines of special high-speed airfoils. 


airfoil. Also, at constant airfoil negative camber, the 
C,,, tends to increase as the airfoil thickness increases. 


The effect of combined negative and positive ‘‘S 
ber on airfoil C,,, has probably never been investigated. 


Call- 


However, some special high-speed airfoils have been 
developed? which have negative camber in the forward 
portion of the chord and positive camber over the aft 
fraction of the chord, as shown for the NACA 66(215)- 
316 airfoil, the best of these five airfoils, in Fig. 7. The 
ratios of +C,,,/—C,,, for these airfoils are all greater 
than 1.0—1.e., 1.22, 1.31, 1.42, 1.32, and 1.21. Al- 
though these data are not conclusive, they suggest that 
positive camber near the trailing edge is more import- 
ant than positive camber near the leading edge for high 
Cy, or, conversely, negative camber near the trailing 
than negative 
camber near the leading edge. it is not likely 
that the “‘W”’ camber airfoils advocated in this paper, 
which have positive camber at about the 85 per cent 
chord station, will have appreciably lower values of low 
speed C;,, than symmetrical airfoils of the same thick- 


edge has a greater adverse effect on C),, 
Hence, 


ness. 

In high-speed airplane design, the trend has been 
toward thinner symmetrical airfoils in order to obtain 
higher M/,’s._ There has been little hesitancy to accept 
the consequent loss in low-speed C,,, and increased 
wing weight resulting from thinner wings in order to 
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obtain the improved transonic aerodynamic character. 
istics. Inverted 
structural thickness and allow a lighter wing than a 
thinner symmetrical airfoil having the same upper sur- 
face curvature, C;,,,and M,. 


-amber airfoils should provide greater 


Some unpublished experimental results indicate that 
“delta’”’ wings the low speed C,,, is in- 


creased by the addition of a sharp leading edge at Rey- 


for sweptback 


nolds Numbers approaching full scale. 

As one of the major causes of the decrease of C,,, re 
sulting from negative camber is the sharper leading 
edge, it appears that this disadvantage of negative 
camber airfoils may not apply to sweptback wings. 
Hence, it is possible that the only argument against the 
use of negative camber for straight wings will ultimately 
prove to be an argument in favor of negative camber 
when applied to sweptback wings. 


Maximum Lift at High Subsonic Mach Numbers 


In the range of M1) = 0.6 to 0.8, the typical positive 
camber airfoil yields a C,,, of about 0.9 to 1.0 for o’s 
of about 8° to 12°. 
data of C,,, of negative camber airfoils at high sub- 
My, but the trends of the negative lift data of 
positive camber airfoils suggest that higher values of 


Unfortunately, there are no test 


sonic 


C),, at higher values of J) may be obtainable with nega- 
tive camber than with positive camber. This possi- 
bility is consistent with the increase in upper surface 
M,, resulting from negative camber, which defers the 
\/y's. This negative 


may be enhanced by upward de- 


compressibility stall to higher 


camber effect on C,,, 
flection of a plain trailing-edge flap, which increases 
the negative camber and further decreases upper sur- 
face velocities. 

The availability of higher C;,,'s at high subsonic / 
resulting from negative camber is highly desirable for 
dive recovery. It should be the subject of immediate 
experimental investigation. 
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NEGATIVE 
A PROPOSED ““W”’ CAMBER AIRFOIL 


Ideally, the condition for maximum J/, is a constant 
upper and lower surface chordwise pressure distribu- 
tion.2 This condition is provided, for an airfoil of zero 
thickness, by the N.A.C.A. a = 1.0 mean line,'! which 
yields a uniform chordwise load distribution. However, 
for an airfoil of finite thickness, the N.A.C.A. a = 1] 
mean line will not have the highest possible J/, for a 
given //¢ and design Ci, because of the pressure distribu- 
tion induced by the thickness form.? 

The maximum J/, for an airfoil of finite thickness 
is obtained for a given C, by the greatest possible chord- 
wise extent of uniform upper surface pressure rather 
than by a uniform chordwise load distribution.” The 
high lower surface positive pressures aft of the stagna- 
tion point contribute to the airfoil lift, but, if uniform 
load distribution is maintained at the leading edge, an 
appreciable increment of negative pressure lift is not 
used at the upper surface leading edge. This lift must 
be provided elsewhere on the airfoil, resulting in higher 
local negative pressures aft on the airfoil and a corre- 
sponding lower J/,. 

To achieve the greatest chordwise extent of uniform 
upper surface pressure, negative camber should be used 
from the leading edge aft to about x/c = 0.60, with 
positive camber aft to the trailing edge. This is the 
essential characteristic of the 
—1.2 


1.5f 


ja = 0.6, Ci 


NACA 66(215)-316 : 
N AC \ 16 ( 15) 3] : la = 1.0, Cy 


I 


airfoil,” which carries no design load on the forward 60 
per cent of the chord and all the design lift aft of x/c = 
0.60, the thickness distribution minimum pressure point. 
The resulting camber line has an “‘S’’ shape, as shown 
in Fig. 7. The sonic MW variation with C;, was esti- 
mated from low-speed test pressure distribution data.” 
The variations are the opposite of those for positive 
camber airfoils, verifying the trend suggested by the 
limited amount of test data for positive camber air- 
foils at negative C,’s. That is, the sonic JJ increases 
with increasing positive C, up to the discontinuity at 
which the leading-edge shock determines the airfoil 
M,. The Ciy’s and C,’s in the low drag range for these 
“S” camber airfoils are approximately the same as for 
the NACA 66 series airfoils having the same thickness 
and approximately the same effective C), with the uni- 
form load, a = 1.0, mean line?. 

Unfortunately, carrying the design lift only over the 
aft portion of the ‘‘S’’ camber airfoils? produces iarge 
negative C,,,’s of the order of —0.050 to —0.100, which 
are prohibitive from the standpoint of longitudinal 
“tuck under” instability. Probably for this reason, 
these airfoils were never tested at high M,’s, so that 
no recommendations concerning their use at high speed 
were made by the authors.” 

On the basis of this information, it is deduced that 
the trailing edge of the camber line should be reflexed 
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TABLE | 
Two Camber Characteristics of N.A.C.A. Mean Lines for Design 
C, = 1.0 
a Maximum y, x/c of Maximum y, Cine/4 
0 6.277 0.30 —0.083 
0.1 6.449 0.30 —(0. O86 
0.2 6.789 0.35 —(0.094 
0.3 7.175 0.35 —0.106 
0.4 7.4389 0.4 —().121 
0.5 7.490 0.45 —().139 
0.6 7.370 0.47 —(0).158 
0.7 7.155 0.50 —().179 
0.8 6.790 0.50 —(). 202 
0.9 6.290 0.50 —() 225 
1.0 5.515 0.50 —(). 250 
with negative camber in order to obtain the necessary 


positive C,,. Hence, a ‘““W’’-shaped camber line 
evolves. 

It is not possible to obtain a higher /, than that of 
the basic thickness form through negative camber 
over the forward portion of the airfoil and obtain a posi- 
tive C,, with a positive C,,. Consider the N.A.C.A. 
mean lines,' three principal parameters of which are 
listed in Table 1 for a C,, = 1.0. As a increases from 
0 to 1.0, the negative C,, increases for positive C;,'s, and, 
as a increases from 0 to 0.5, the maximum camber in- 
But, as a increases from 0.5 to 1.0, the maxi- 


To obtain a positive C,, and 


creases. 
mum camber decreases. 
a positive C;,, an a = 1.0 camber line with negative 
C,, could be used in combination with a camber line 
of a less than 1.0 but having a slightly larger positive 
C;,. The combination of two such camber lines pro- 
duces positive camber in the forward portion of the air- 
foil, because the camber for a less than 1.0 is greater 
than for a = 1.0 and a larger positive C,,; was required 
for a less than 1.0 than for a = 1.0 in order to achieve 
a positive net C;,. Hence, the .1/,’s would be less at 
high-speed C;,’s than those for the basic thickness form 
because of the increased upper surface curvature. 

The only one of the three airfoil design parameters 
positive C,,, positive C,,, and negative camber—which 
is not essential to the airfoil’s aerodynamic high-speed 
characteristics is the net positive C,,. Hence, it is 
advocated that transonic airfoils be developed with net 
negative C,,'s, ““‘W’’ camber, and a positive effective 
“design lift coefficient’ obtained by a positive AC, 
resulting from angle of attack. The AC, provides 
the leading-edge negative P necessary to obtain the 
greatest possible chordwise extent of minimum uni- 
form P. The negative camber in the forward portion 
of the chord relieves the usual mid-chord maximum 
negative P at high-speed C,’s. The positive camber at 
about x/c = 0.85, the center of the ‘““W’ 
extends the upper surface negative P aft as far as pos 
sible. Finally, the trailing-edge negative camber, 
the second half of the “‘W”’ camber line, provides a small 
trailing-edge down load sufficient to yield a net posi 


camber line, 


tive Cp. 
As previously discussed, the leading-edge peak nega 
tive P causes a shock of low intensity for a considerable 
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Fic. 9. Lift and drag divergence Mach Numbers of NACA 
847B110 airfoil at positive and negative lift coefficients. 
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Section drag polars for NACA 847B110 airfoil at 
several Mach Numbers. 


Fic. 10. 





Section pitching moment curves for the NACA 847B110 
airfoil at several Mach Numbers. 


range of C,, which does not determine the airfoil diver- 
gence M, until it extends over about 20 per cent of the 
chord at higher C,’s. 


The Detail Design of ‘‘W"’ Camber Transonic Airfoils 


The ““W"’ camber concept has been applied to several 
airfoils derived with various combinations of N.A.C.A. 
thickness forms and camber lines. The best two of 
these airfoils are presented below. 

All the airfoils were arbitrarily developed with a net 
C,, = —O0.1, and then a AC, = +0.2 was added to 


i 
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yield a net working C,; = +0.1 at My) = 0. Also, all 
airfoils were calculated with maximum ¢t/c = 0.10. 


Large theoretical positive C,,'s were maintained, be. 
cause test experience indicates that the theoretical 
trailing-edge load is not realized physically. 

One of the two airfoils is designated in the N.A.C.A. 
system as the 


(2 = 0.7, C, = — 0.3) 
NACA 67(015)-(-1)10 <@ = 0.8, C, = +0.6 
(1k C.~ -O8 

which has a theoretical net C,, = +0.0325 at the net 


C,, = —0.1. The theoretical pressure distribution at a 
working C, = +0.1, the airfoil contour, the tabulated 
airfoil coordinates, the camber line ordinates and slope, 
and velocity and pressure distributions are shown in 
Fig. 12. 

The V/, for first occurrence of the pressure distribu- 
tion shock discontinuity, estimated from correlated 
flight data,* is noted as 0.864 on the lower surface, 
At the 7, = 0.864, the C; at constant a is estimated 
to be 0.20 instead of the C,; = 0.10 at My = 0. The 
estimated 7, was based on the maximum negative 
P atx/c = 0.70 rather than on the higher negative P 
at the leading edge. 

The low-speed maximum negative P = —0.253 on 
the lower surface is estimated*® to become shocked at 
My = 0.864, while a low-speed P = 0.15 yields local 
sonic velocity® at JJ = 0.864, showing that the local 
flow is slightly supersonic before a shock discontinuity 
obtains in the pressure distribution (Fig. 4). 

The characteristics of a second airfoil, designated 


“<0, C= -02 
' o=0@8 C, = +62 
NACA 16-(- 
ee I ee Cw ke 
o—ik C= 08 


are presented in Fig. 13. The net theoretical C,, at 
C, = —0.1 and My = 0 is +0.0233. 
\/.for first shock is 0.865. 

The two airfoils presented herein are considered 


The estimated 


representative of the design philosophy and are not 
necessarily the optimum configuration. It is probable 
that a better combination of negative C),, chordwise 
camber distribution, and AC,, due to a can be found 
to yield an optimum transonic airfoil. 

The writer has not been able to obtain high-speed 
test data on the above airfoils, but there are unclassi- 
fied test data available which can be presented as ex- 
perimental proof of the design principle. The specially 
developed NACA 8 series high-speed airfoils’ possess a 
camber line having the shape of the letter ““M.”’ In- 
verting the 8 series airfoils yields an airfoil with the ad- 
vocated ‘““‘W”’ camber distribution. The best of the 
NACA 8 series airfoils is the NACA 847B110, the cam- 
ber line of which is shown in Fig. 7. The lift and drag 
divergence 1/,’s for this airfoil are shown in Fig. 9 as a 
function of positive and negative C,’s with the nega- 
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NEGATIVE 
TABLE 2 


Per Cent 
Inferiority 


M (Cy for Cu, ., Ca, in Ca at Equal + C of + C 
0.80 —0.11 0.0071 0.0098 38 
() 82 —().19 0.0074 0.0182 146 
0.85 ~0.16 0.0202 0.0338 65 


tive C, range mirrored in superposition on the positive 
C, range to accentuate the visual comparison of positive 
and negative camber. It will be seen that the .1/,’s in 
the negative C; range are considerably higher than those 
in the positive C, range. 

As the true criterion of negative camber drag super- 
jority is the absolute drag, the section drag, Cy, has 
been plotted for the NACA 847B110 airfoil in Fig. 10 
as a function of positive and negative C, at MJ) = 0.3, 
0.75, 0.80, 0.82, and 0.85. 
as My increases, the C, for minimum C, becomes nega- 


These curves show that, 


tive. For the three ‘higher M's above, a comparison 
of the minimum C, at a negative C,, with the corre- 
sponding Cy at the equal positive C), is presented in 
Table 2. 

The variations of C,, with C, at constant ./) for the 
majority of the 8 series airfoils’ indicate irregularities 
and double inflections, but the NACA 847B110 airfoil 
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Characteristics of the 
42 = ().7, Cc, = —0.3 
NACA 67(015)-(-1)104 a = 0.8, Ci, = +0.6 2 airfoil 
la =10, CG, = 0.44 
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40.002 | 4364 | 39998 |5376 |-4960/-0005 |: 119 |-25: | 1095 |-200 
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Fic. 13. Characteristics of the 
a = ().7, C = 0.3 
“AC » 220.8, C, = +0.2{ .,. 
NACA 16-(-1)105 ‘ ; : “ air 
4 (oon aan iirfoil 
a=10, G = —0.4 


yields quite regular and continuous variations of C,, 
with C, and J/y. It is interesting to note, in Fig. 11, 
that, as JJ, increases, the C,,, approaches zero nega- 
tively and then equals zero at JJ = 0.80 and 0.85. 
Also, the zero lift line is coincident with the chord line, 
indicating that the positive and negative camber are 
mutually compensating.° 

Another specially developed high-speed airfoil that 
incorporates unconventional chordwise camber dis- 


tribution is the 


I| 


en ey (a= 1.0, C,, 0.6) 
aNd d 66 09)-2 Z . 
emda ee 


airfoil.* * A comparison of this camber line with the 
a = 1.0 camber line of the NACA 66-210 airfoil is made 
in Fig. 7. The greatly reduced camber in the forward 
portion of the NACA 66 9)-210 airfoil compared to the 
NACA 66-210 airfoil reduces the upper surface curva- 
ture and, consequently, increases the ./,. This effect 
is supplemented by the high loading on the trailing-edge 
portion aft of the thickness distribution minimum pres- 
sure point. 

A comparison of the C, and C, divergence \V/,’s for 
these two airfoils is shown in Fig. 8. It can be seen that 
the divergence /,'s for the NACA 66 09)-210 airfoil are 
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slightly higher than those for the NACA 66-210 airfoil 
in the positive C; range, but this increase is obtained at 
the expense of a large increment in negative C,,,, which is 
intolerable because of the ‘‘tuck under’’ longitudinal 
instability induced thereby. 

It should be noted that a greater increase in diver- 
gence M/,’s over an extremely adequate range of C;'s 
for high-speed flight is obtained by merely inverting 
the basic NACA 66-210 airfoil than is obtained by using 
the special NACA 66,09)-210 airfoil. 

Some significant data* are reproduced in Fig. 8 
which contrast the much greater range of C, for 
C,and C, high divergence M, realized in test, with the 
lesser range of C, for high A/, predicted from low-speed 
pressure data or theoretical pressure distributions. 
Also, the Cy divergence /,’s are higher than the /,’s 
predicted from low-speed test and theoretical pressure 
distribution data. These effects are believed to be 
produced by the following three effects: 

(1). The local A/ for shock on an airfoil is slightly 
supersonic instead of sonic, as assumed in the von Kar- 
man relation.® 

(2) The low-intensity premature leading-edge shock 
is absorbed by the chordwise aft subsonic flow, and has 
no effect on the airfoil C,, Cz, and C,, characteristics, at 
small a’s. 

(3) The test C, and Cy divergence M,’s are plotted 
against test C, at the divergence J/,, while the theoreti- 
cal M,’s are plotted against the low-speed C; unin- 


creased by compressibility. 


CONCLUSIONS 


(1) Airfoil divergence \/,’s in the high-speed C, 
range can be increased by the use of negative camber. 

(2) The maximum J/, is obtained with negative 
camber over the forward chordwise portion of the air- 
foil and positive camber aft to the trailing edge but at 
the expense of large negative C,,’s. 

(3) The best compromise camber distribution has 
the form of a ‘“‘W.’’ The negative camber in the for- 
ward portion of the chord relieves the mid-chord maxi- 
mum negative P at low a. The positive camber at 
about x«/c = 0.85, the center of the “‘W,”’ extends the 
upper surface negative P as far aft as practical. The 
trailing-edge negative camber, the second half of the 
“W,”’ provides a small trailing-edge down load sufficient 
to yield a net positive C,,,,. 

(4) The NACA 847B110 airfoil, when inverted, has 
the proposed ““‘W”’ camber distribution. The nega- 
tive lift transonic test data for this airfoil provide test 
verification of the advantages of the proposed “‘W”’ 
camber line. 

(5) Negative camber airfoils should provide greater 
structural thickness than thinner symmetrical airfoils 
having the same upper surface curvature C,,, and diver- 


gence /,’s. 





1951 
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(6) Negative airfoil camber provides the positive 
Cy, required to trim every stable airplane without con- 
trol deflection at positive C;’s. 

(7) Negative camber should tend to eliminate the 
transonic tendency toward longitudinal instability 
known as “‘tuck under” caused by dC,,/dC, and Cy, 
variation with Mo. With negative camber, the C,,, 
variation tends to compensate for the dCy/dC,, vari- 
ation. 

(8) Two “‘W” camber airfoils presented are desig- 
nated in the N.A.C.A. system as the 


a=0.7, C, = —03 
NACA 67(015)-(-1)10 <a = 0.8, C, = +0.6 
a=1.0, C, = —0.4 
and the 
a=0.7, C, = —0.3} 
9 


NACA 16-(-1)10 Ja = 0.8, C, = +0.2 
a=0.9, C, = +0.4 
a=1.0, C, = —0.4 


airfoils. At AJ) = O, the theoretical values of C,, at the 
net C;, = —0.1 are +0.0325 and +0.0233, respectively. 
The estimated J7/,’s for first occurrence of the pressure 
distribution shock discontinuity on the lower surface 
are 0.864 and 0.865, respectively, at the corresponding 


C; i +0.2. 
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High-Efficiency Supersonic Diffusers’ 


E. P. NEUMANN? anpb F. LUSTWERKt!? 
Massachusetts Institute of Technology 


ABSTRACT 


Supersonic diffusers having a variable throat were designed 
with wind-tunnel application in mind. The test results indicate 
that a stagnation pressure recovery of 79 per cent at AJ = 2.22 
and 53 per cent at J = 2.92 can be obtained without boundary- 
layer suction. These values can result in an operating power 
reduction of 60 per cent over conventional design at WJ = 2.22. 

The effect of scale on diffuser performance is discussed. De- 


sign parameters for a diffuser design are given. 


NOMENCLATURE 


a@ = cross sectional area, sq.ft. 

D = diameter of test section 

D, = equivalent diameter; 4 (cross-sectional area perim- 
eter 

g = acceleration given to a unit mass by unit force, ft. per sec. 
absolute 

k = ratio of specific heats, for air k = 1.40 

L = length of test section 

M = Mach Number = V/V gkRT 


p = pressure, lbs. per sq.in. absolute 


R = gas constant (1,545 ft./lb./°F. Ib. mole divided by 
molecular weight 

T = temperature, °F. absolute 

V = mean velocity of fluid stream at given cross section, ft. 
per a 

8 = oblique shock angle 

6 = wedge angle 

» = diffuser efficiency, work of isentropic and adiabatic 
compression between initial condition and final pres 
sure divided by kinetic energy expended [see Eq. 
{ 

Subscripts 


a refers to state of fluid stream after first oblique shock (see 
Fig. 7 

b_ refers to state of fluid stream after second oblique shock (see 
Fig. 7) 

c refers to minimum cross section of diffuser 

section at minimum area of accelerating 


t refers to cross 


nozzle 

0 refers to stagnation state of fluid stream at entrance of 
accelerating nozzle 

1 refers to state of fluid stream at the diffuser entrance 

2 refers to state of fluid stream downstream of ‘‘transverse”’ 
shock where stream is subsonic 

3 refers to state of fluid stream after diffuser where the ve- 

locity is approximately zero 

Received September 1, 1950. 

* This work was carried out as part of the Guided Missiles 
Program at the Massachusetts Institute of Technology and was 
sponsored by the Navy Bureau of Ordnance. 
Figs. 2, 10, and 11 were obtained from Krenkel.* 
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INTRODUCTION 


ip EFFICIENT DECELERATION of a supersonic air 
stream in the presence of boundary layer is an im- 
portant problem in the design of wind tunnels, com- 
pressors, intakes to high-speed aircraft, and other similar 
devices. For this reason an investigation of this prob- 
lem was undertaken at the Gas Turbine Laboratory at 
M.1.T. several years ago. The first results of this work 
were reported by the authors! in a paper that discussed 
efficiency and operating characteristics for diffusers of 
fixed geometry. These diffusers were not designed to 
take advantage of the theoretical gain that might be 
obtained by decreasing the throat area of the diffuser 
after the shock had passed through. In the present 
paper, the earlier work is extended for the case where 
the diffuser walls are made adjustable to permit reduc- 
tion of the throat area after starting—that is, after the 
shock has passed through the diffuser throat. 

In Meteor Report No. 13,' separation of the stream 
from the passage wall was observed to accompany a 
“normal” shock. However, the change in state across 
the separated region was approximately the same 
as that predicted by a 
analysis, provided the separated region was contained 


simple one-dimensional 


in a constant-area The analysis assumed 
that wall forces in the direction of flow within the con 
That this assump- 


passage. 


stant-area passage were negligible. 
tion agrees adequately with reality was confirmed. The 
optimum length of this constant-area passage was de- 
termined experimentally to be 8 to 12 equivalent 
diameters for a range of Mach Numbers from 1.8 to 
4.2. The results of Meteor Report No. 13 are included 
with the new results in order to facilitate comparison. 
The work reported in the present paper on variable 
geometry diffusers started with the simple concept that 
the shock and attendant separation could be adequately 
handled from a pressure recovery standpoint and that a 


subsonic diffuser of high efliciency could be employed 


— — 
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Fic. 2. Apparatus for central-body diffuser tests. 
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Apparatus for movable sidewall diffuser. 


downstream of the shock. With this concept it ap- 
peared reasonable to expect high overall efficiencies 
for a variable-geometry supersonic diffuser such as 
might be employed in a wind tunnel. After several un- 
successful attempts to achieve a variable-diffuser throat 
and high-efficiency performance? by employing designs 
similar to those shown in Fig. 2,* the diffuser shown in 
Fig. 3 was designed and tested. With this diffuser a 
stagnation pressure recovery of 79 per cent (diffuser 
efficiency, 87 per cent) was obtained at a Mach Number 
of 2.22. 


* Figs. 2, 10, and 11 were obtained from reference 2. 
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One might consider what this adjustable diffuser 
could mean in terms of the power required to operate 
a wind tunnel at the Mach Number mentioned. For 
comparison, the conventional diffuser usually employed 
in wind tunnels is shown in Fig. 5, described as a Type 
I diffuser. A reduction in power requirement of 60 
per cent could be achieved in the limit with the adjust- 
able diffuser. If one thought in terms of a large tun- 
nel requiring 100,000 kw., such a saving would be signifi- 
cant. 


EXPERIMENTAL APPARATUS 


In an attempt to expedite the collection of data and 
minimize the effort required to construct experimental 
apparatus, small-scale equipment has been used 
throughout. 

Fig. | shows the arrangement of the test apparatus. 
High-pressure air, 125 Ibs. per sq.in. gage, at a rate of 
0.55 Ib. per sec. was available from the air compressor. 
The apparatus used to provide movable walls in the 
Two small-scale supersonic- 


) 


diffuser is shown in Fig. 3. 
wind-tunnel nozzles were used to provide the supersonic 
stream at the diffuser entrance. Fig. 4 shows one of 


these nozzles. 


COMPARISON OF DIFFUSER TYPES 


Experimental results are summarized in Figs. 5 and 
6. In order to facilitate comparison of the variable- 
geometry diffuser with some other types of diffusers, 
the test results previously obtained by the authors! 
are presented here. 

Fig. 5 shows measured values of the percentage of 
stagnation pressure recovered plotted vs. inlet Mach 
Number. The percentage of stagnation pressure 
recovered is determined by measurement of initial 
stagnation pressure and temperature (/») and (7%) (see 
Fig. 7), final static pressure (p3), the static pressure at 
the diffuser entrance (f;), and the weight flow (w) 
through the diffuser. The weight flow is computed 
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Po 2g ( k )( 2 we =} 
OU VT, UR \b + 1/\R + ) 


w = 0.531a,(po/V T») (1) 
The per cent stagnation pressure recovery 


, (k — 1) 
= ps (100)/pi: | M2? ———— + 1 


k/(k 1) 


(100) 
Istag 
(2) 
The value of the Mach Number at the diffuser entrance 
(M,) is determined by 


—1+ 


bo\*/ai\7/k — 1 9 2/(k — 1) 
Le ( ( ) ( )( (3) 
PiJ \ai/ \k + 1/\R + 1 


where a, and a; are measured values of the accelerating 
nozzle throat area and the diffuser entrance area, respec- 
tively. 

An alternate method for presenting these data is 


shown in Fig. 6. Here, the diffuser efficiency (7) is de- 
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fined as the ratio of isentropic increase in enthalpy to 
the change in kinetic energy in the diffuser. The isen- 
tropic increase in enthalpy is computed between the 
state preceding the diffuser entrance and the static 
pressure after the exit to the subsonic diffuser. The 
efficiency may be expressed in terms of symbols as fol- 


lows: 


_ [2 (k — 1)] [(ps/p.)“ D/e 1] 
" M,? — (Vs?/gkRT:) 4) 


where V3 is the velocity leaving the subsonic diffuser. 
For the tests reported, V; was kept small. Had the 
diffuser been charged with this leaving kinetic energy, 
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Fic. 9. Summary of analytical results. 
the efficiencies reported in Fig. 6 would in all instances 
be within 1.5 per cent of those shown. 


ANALYTICAL RESULTS 


The previous work! indicated that a ‘normal’ shock 
occurs during the starting process. The minimum 
starting area, as computed from a one-dimensional anal- 
ysis that neglects friction, agrees with the actual 
data within the limits of experimental accuracy. 


By decreasing the area of the diffuser throat after 
the “‘normal’’ shock has passed through, the power re- 
quired for the operation of a wind tunnel can be de- 
creased. Without stability considerations, the Mach 
Number could be decreased to 1. 
and friction considerations limit the minimum Mach 


However, stability 


Number to a larger value. 


Theoretically, it is possible to achieve a reversible 
compression to the sonic velocity by using an infinite 
number of infinitesimal shock waves, as proposed by 


Oswatitsch. Practically, it is necessary to take into 


SCIENCES—JUNE, 1951 
account the irreversibilities due to finite oblique shocks 
and the limitations imposed by the boundary layer. 

For the purpose of this report, only wedge-shaped 
entrances are considered. The model selected for anal- 
ysis is shown in Fig. 8. It is assumed that the dif- 
fuser throat has been reduced to area a, after the nor- 
mal shock has passed by and that frictional forces are 
negligible. The value of a, is defined by the inter- 
section of the oblique shock with the wedge of angle 
6 at point A. 
actly canceled by the expansion wave at point A. A 


That is, the incident shock will be ex- 


normal shock is then assumed to occur at the area a, 
Following this shock, the stream is decelerated to zero 
velocity in a reversible subsonic diffuser. 

Using this model, one can compute the pressure ratio 
(p3/pi) for a given entrance Mach Number and wedge 
angle @ from conventional equations.’ The results of 
such computations are shown in Fig. 9 and are also 
summarized in Fig. 6. 

For the experimental model, a wedge angle @ of 10 
was chosen. It can be seen, however, from Fig. 9 that 
this angle gives a diffuser efficiency of 93.8 per cent, 
which is 1.6 per cent less than the minimum efficiency 
predicted at 6 = 12°. The smaller angle with the 
slightly lower efficiency was selected because of sta- 
bility considerations. Note that JJ, = 1.48 at @ = 
10°, whereas M7, = 1.32 at @ = 12 When the Mach 
Number in the throat is near 1, a small disturbance 
propagated upstream through the subsonic diffuser can 
cause complete flow breakdown—that is, the ‘‘normal” 
shock will move upstream through the converging 
section to the diffuser entrance. Should breakdown 
occur, the area a, would have to be increased until start- 
In addition to 


ing conditions were again established. 





Fic. 10. Schlieren photograph of central-body diffuser. Flow 
is from left to right. The shock from the 5° half-angle wedge is 
canceled at point A by the expansion around the 6° corner. 
A weak expansion wave leaves the corner. Effects farther down 
stream are caused by the back pressure so that the ‘‘normal” 
shock region is seen. 





Fic. 11. 


Schlieren photograph for conditions similar to those 
of Fig. 10 except that the central wedge has been moved down- 
stream so that the shock and expansion wave no longer com- 
pensate. 
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HIGH-EFFICIENCY 


Cross-Section 


Dimensions D. 
In.) (In.) M 
18 x 24 20.6 2.0 
9x8 8.47 1.9 
1.292 XK 1.292 1.292 2.22 
l-in. diameter l 2.05 
0.633 X 0.633 0.633 2.22 
l-in. diameter l 2.2 


stability considerations, it is desirable to select an angle 
@ smaller than the optimum to make allowance for in- 
Note that the effi- 
ciency of a “‘normal’’ shock at 1/; = 1.5 is 92 per cent, 


efficiency in the subsonic diffuser. 


whereas a subsonic diffuser of 90 per cent efficiency is 
considered good. Consequently, frictional effects in 
the subsonic diffuser, not considered analytically, will 
cause the best diffuser efficiency to occur at a somewhat 
lower value of 6 thart is indicated from Fig. 9. 


EXPERIMENTAL RESULTS 


The analytical model described above was selected 
as a basis for comparing measurement and theory. 
Some of these data are shown in Fig. 9. It is interest- 
ing to note that the measured value of diffuser efficiency 
is 94 per cent of the computed value at JJ, = 2.22 and 


90 per cent at JJ, = This comparison is based 


on a 10° wedge angle 6. A diffuser section was used 
which was similar to the Type V diffuser (Fig. 5), ex- 
cept that the angle in the “normal” shock region was 
1.2° instead of 50 min. The 1.2 


in an effort to improve starting by allowing for bound- 


divergence was made 


ary layer. 

‘This analysis assumes that the incident shock wave 
is canceled by the expansion wave at point A. Fig. 
10 shows experimental conditions comparable to this 
The model used to obtain this picture 
half-wedge angle, 


assumption. 
employs a central wedge with a 5 
whereas the angle between the stream and the wall is 
6°. With this choice of angles the shock wave inci- 
dent at the corner is weaker than the expansion wave; 
consequently, it will be canceled, and a weak expansion 


wave, equivalent to a 1° turning angle, will leave the 


corner. 
One can see in the photograph that a weak 
expansion wave leaves the corner. The boundary 


layer at the wall is turbulent. Fig. 11 shows condi- 
tions similar to those of Fig. 10, except that the central 
wedge has been moved downstream so that the shock 
and expansion wave no longer compensate. It can be 
11 that the interaction of shock wave 
layer is an important consideration in 
diffuser design. Note in particular the separation phe- 
nomena in the region where the reflected shock is 
(See reference 14 for a 


seen from Fig. 
and boundary 


incident on the central body. 
discussion of shock-wave and boundary-layer inter- 


action. ) 


SUPERSONIC 


TABLE 1 


DIFFUSERS 
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Mo Corrected 
( Approx.) L/D, to M =2 
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Fic. 12. Multicellular diffuser. This type of construction was 


employed to reduce the overall length of the diffuser. 


SCALE EFFECTS 


The majority of the experimental data reported here 


are for small-scale tests 
sections around 1!/, by 1'/, in. and circular cross sec- 


It is interesting to note that 


namely, rectangular cross 


tions of l-in. diameter. 
the same results were obtained with the rectangular and 
circular sections, provided data were reduced to equiva- 
lent hydraulic diameter. 

Additional data have been obtained on two diffusers 
employed in larger supersonic wind tunnels at M.I.T. 
with test sections of 9 by 8 in. and 18 by 24 in. The 
diffusers used were of the Type II, shown in Fig. 5. 
The results of these tests are shown symbolically as 
large solid squares on Fig. 5. These results indicate 
that the small-scale results can be extrapolated readily 
to larger scale. The only scale effect noted is that the 
required length of the constant-area section seems to 
be somewhat shorter (three diameters shorter at JJ = 
2 for the 18- by 24-in. tunnel). A summary of the 
data on scale effect is shown in Table 1. The length of 
shock region seems to increase as Reynolds Number 
VD.p/ mo decreases. 


LENGTH REQUIREMENTS FOR SUPERSONIC DIFFUSER 


The length requirement for a supersonic diffuser is 
often an important consideration from the standpoint 
of application—for example, building space in a wind- 
tunnel installation and length of aircraft intakes or 
passages in supersonic machinery. With these con- 
siderations in mind, some experimental data were ob- 
tained on a diffuser of the type shown in Fig. 12. 

The required length for the constant-area section was 
about nine and one-half hydraulic diameters for an 
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entrance Mach Number of 2.22. This observation 
bears out other data on scale effect—namely, the re- 
quired length of the shock region for a given Mach 
Number depends primarily on the equivalent length- 
diameter ratio (L/D,). 


A significant result obtained from the tests on the 
diffuser of the type shown in Fig. 12 is that the overall 
length can be reduced by 25 per cent over the Type II 
diffuser of Fig. 5 without loss in efficiency. The meas- 
ured stagnation pressure recovery was 61 per cent and 
60 per cent, respectively. Further subdivision of the 
shock region might further reduce the length require- 
ments. 


CONCLUSIONS 


A small-scale variable-throat diffuser was designed 
for wind-tunnel application which has a stagnation 
pressure recovery of 79 per cent at M = 2.22 and 53 per 
cent at MW = 2.92. Boundary-layer suction was not 
necessary to achieve these values. The starting charac- 
teristics of this diffuser were as good as any other type of 
fixed-throat diffuser tested. On the basis of tests to de- 
termine the scale effect on fixed-geometry diffusers, it 
appears reasonable to assume that the variable geom- 
etry results can also be extrapolated to larger scale 
application. 
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Dynamic Analysis of Aeroelastic Aircraft 
by the Transfer Function-Fourier Method’ 


JAMES B. REAt 
Douglas Aircraft Company, Inc. 


SUMMARY 


The purpose of the paper is to outline a method for determin- 
ing in a single analysis the aircraft flight-path stability, the 
dynamic tail loads, and the flutter response of the system. This 
method treats the airplane itself as a servomechanism and is 
called the Transfer Function-Fourier (TFF) method. 

The TFF method is first explained in general, followed by 
simple examples to demonstrate the concept of the Transfer 
Function and the use of the Fourier technique in obtaining the 
transient response. AISo included is a brief summary of the 
history leading to the development of the TFF method in the 
field of automatic control. Examples are then given to show the 
conventional application of this method to the dynamic stability 
analysis of simple feedback servomechanisms and automatic 
control systems. 

The paper then points out that, through extension of these 
fundamental concepts, an aeroelastic aircraft may be analyzed 
dynamically by treating it as a multiple loop servomechanism 
with aerodynamic and aeroelastic feedbacks. This allows many 
of the methods and “‘tools’’ that were developed years ago by the 
electrical-servomechanism engineers in their study of the dy- 
namic stability of communication networks, feedback amplifiers, 
servomechanisms, and automatic control systems to be applied 
directly by the TFF method to the dynamic analysis of aero- 
elastic aircraft. The paper also points out that many of the an- 
alysis techniques recently developed in the study of the dynamics 
of aeroelastic aircraft by conventional methods, especially the 
matrix techniques developed in the field of flutter, may be applied 
to the dynamic analysis of systems for automatic control. 

Practical examples are then given to demonstrate the ex- 
tension of the TFF method to the dynamic analysis of rigid-body 
and aeroelastic aircraft. These examples include: 

(a) A rigid-body longitudinal dynamic analysis of an airplane 
to determine the transient maneuvering tail loads developed by 
elevator displacements. 

(b) A longitudinal dynamic stability and performance analy- 
sis (and a limited flutter analysis) of an airplane equipped with 
a mechanical gust load alleviation system that is actuated by 
elastic wing deflection. 

(c) A longitudinal dynamic analysis of an airplane to deter- 
mine by a single analysis: (1) the stability of the aircraft flight 
path, (2) the dynamic tail loads, and (3) the flutter response of 
the control system and the multielastic structure, where these 
responses result from arbitrary gust disturbances or control force 


inputs. (This more complete analysis includes the degrees of 


Presented at the Aeroelasticity and Structures Session, Annual 
Summer Meeting, I.A.S., Los Angeles, July 12-14, 1950. 
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freedom of the airplane as a whole, its control system, its multi- 
elastic structure, and aerodynamic circulation.) 

In conclusion, the paper points out the main advantages of the 
TFF method and makes recommendations for further study and 
development of the method. 


INTRODUCTION 


yy THE AIRCRAFT INDUSTRY TODAY, many separate 
dynamic analyses pertaining to the air frame are 
made. For example, there are: (a) dynamic stability 
analyses of the aircraft flight path resulting from gust 
disturbances or control force inputs; (b) dynamic 
analyses to determine the transient sttuctural loads 
(especially tail loads) resulting from gust disturbances 
or intentional maneuvers; and (c) flutter analyses of 
the aircraft control system and its multielastic struc- 
ture. 

It is the purpose of this paper to outline a method for 
accomplishing these analyses in a single analysis. This 
method treats the airplane itself as a servomechanism 
and is called the Transfer Function-Fourier (TFF) 
method. 

Before demonstrating the application of the TFF 
method to the dynamic analysis of aeroelastic aircraft, 
the method will be explained in detail, showing first its 
application to a simple mechanical system, and then 
its conventional application to simple feedback servo- 
mechanisms and automatic control systems. 


(II) EXPLANATION OF THE TFF MetTHop 


The TFF method can be explained as follows: As- 
sume that the transient response of a linear system is 
required for the case when the input to this system is 
arbitrary inform. The TFF method, as applied to this 
problem, would then consist of (1) computing the 
transfer function of the system (which is nothing more 
than the forced steady-state output-to-input vector 
response ratio for harmonic input); (2) representing 
the actual input by its Fourier equivalent; (3) operat- 
ing on each term in the input by the corresponding 
value of the transfer function, to obtain the Fourier 
equivalent of the transient output; and (4) summing the 
terms in the resulting output to get the net transient 
response of the system to the arbitrary input. 

When the system is of the multiple open loop or the 
multiple closed loop feedback type, it is convenient in 
computing the transfer function of the system, step 
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Fic. 1. 


Simple mechanical system. 


(1), to construct a block diagram and treat the system 
as a servomechanism. 

The TFF method can be further clarified by a more 
complete explanation of the concept of the transfer 
function. 


(a) Transfer Function Concept 


In general, the transfer function of a dynamic system 
is that function that “‘transfers’’ the system input to the 
output—i.e., the output is obtained by operating on 
the input with the transfer function. 

In particular, the transfer function for a linear system 
with harmonic input is of the vector form (A R)o,e" PA); 
where (A R)o; is the amplitude ratio of output to input, 
a function of frequency and (PA)po; is the phase angle 
between output and input, also a function of frequency. 
The linear harmonic transfer function is usually repre- 
sented by plots of (A R)o,; vs. frequency and (PA )o; vs. 
frequency, or by a single polar plot of (AR)o, vs. 
(PA )o, with frequency as a parameter. 

The general transfer function concept is also valid for 
a nonlinear system, except, of course, it does not have 
the form (A R)oye“’4’o;;_ nevertheless, it usually can be 
represented by an output-to-input ratio. 

The mathematical form of the linear harmonic trans- 
fer function can be demonstrated with the usual “‘stock”’ 
example of a simple mass-dash pot-spring system (Fig. 
1) whose physical characteristics can be described by a 
linear differential equation with constant coefficients: 


F = force input 

x = displacement output 

K = spring constant 

C = viscous damping constant 
M = mass of cart 


The differential equation of motion for the system is 
M(d*x/dt®) + C(dx/dt) + Kx = F 
Then, for harmonic force input, 
F = Fe“? 
and the output-to-input ratio, or transfer function, is 


% ] 


Pie. == ; : 
F K(1 — B? + 2¢Bi) 


where 
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F ws forced frequency 
Wn natural undamped frequency 
: C ¢ : ‘ 
5; =f =- - = damping ratio 
Cou. 2V7MK 
1 = /—1 
Thus, 
Ka l 
(A R) = _ = f 
Fo KV (1 — B?)? + (2¢8)? 
(PA) rp = —tan—! [2¢8/(1 B?) | 


where, by definition, 

(TF). = (AR), » es 
In other words, when the input, F, is multiplied by the 
transfer function, (7F),-, the output, x, 1s obtained, as 
follows: 
(Input) (Transfer Function) 


(F)(TF) ip = (Fyei®!) (AR) gp ef? 4) 


i|w;t+(PA) pl 


(Output) = 
x = 
= Xe 


Thus, it can be seen that the transfer function is es 
sentially a vector operator that ‘“‘transfers’’ the ampli- 
tude of the input vector to the amplitude of the output 
vector and shifts the phase angle of the input vector to 
the phase angle of the output vector. 

Once the transfer function for this simple system, 
(TF),», has been determined, it may be combined with 
functions to determine the overall 


other transfer 


transfer function of a larger system [as shown in 
Parts (III) and (IV) of this paper], or it may be used 
to compute the transient response of the displacement, 
x, to an arbitrary input force, F. This is done by 
representing the input by its Fourier equivalent* and 
then operating on each term in the input with the cor 
responding value of the transfer function to obtain the 
Fourier equivalent of the output, which may then be 
summed to give the actual transient response. 

This process of determining the transient response by 
the Fourier technique will now be explained in further 
detail. 


(b) Transient Response by Fourier Technique 


In general, the Fourier technique can be applied to 
two basic problems that are commonly found in the 
application of the TFF method. These problems are: 
given the steady-state frequency response data (i.e., 
the transfer function), find the transient response to a 
given arbitrary input; or, given the input and the cor- 
responding output, find the transfer function. Both 


of these problems have practical application in the 


* Most practical inputs may be represented by a Fourier series 
however, the usual Fourier precautions relating to boundary 
conditions and discontinuities should be observed in applying 
this technique. The most common cases where these precautions 
should be considered are when the transfer function goes to in 
finity at zero frequency or to a finite value at infinite frequency 
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In the 


first case, it is only necessary to represent the input as 


analysis and synthesis of a dynamic system. 


a Fourier series, multiply the amplitude and shift the 
phase of each component as indicated by the value of 
the transfer function for that frequency, and sum the 
resulting components. In the second case, it is possible 
to express both the input and the output as Fourier 
series with the same fundamental frequency, and then 
the ratio of terms of corresponding frequency is the value 
of the transfer function at that frequency. 

It is the first case that is most used in the application 
of the TFF method; however, since both cases are im- 
portant, the theory covering both will be given. The 
theory for the first case is as follows: 

Any ordinary periodic function can be represented by 
a Fourier series. If the independent variable is ¢ and 
the function is f(t), then the exponential form of the 


series is 


n + x 
fjy= > Cine (1) 

n=-—« 

where 

| my 4D 

C(n) = +f f(t) ee" dt (2) 
lo 

forn = 0, +1, =2, =3,..., and w = 2x/T, where 


T is the period of the function f(t) and f) is the initial 
value of f. 

The exponential form of the series is used here only 
for convenience in the following development. 

Nonperiodic functions can be handled by use of the 
Fourier integral or the Laplace transform, but it will be 
sufficient for the present to consider only a part of the 
nonperiodic function. If the part is long enough for the 
function to approach closely its steady-state value, it 
will contain all the pertinent information. This useful 
part of the function can be expanded in a Fourier series 
in a number of ways; the following is convenient and 
has a simple physical interpretation: Take the funda- 
mental period to be 7, such that the transient has 
damped to a satisfactory small amplitude at a time less 
than (1/2)7. To complete the cycle, invert the first 
part and shift it vertically and to the right so that it 
joins the first partial cycle. The process is illustrated 
in Fig. 2 and means physically that the moving element 
is imagined to return to its original position after the 
transient has disappeared, and in the manner in which 
it was originally displaced. If the period T is divided 
into two equal parts, half of the coefficients will be non- 
zero because of the symmetry. To get both odd and 
even harmonics, the two parts of the cycle would have 
to be taken of different duration. 


Definitions 
I(t) = input to system 
O(t) = output from system 
(AR)o, = amplitude ratio of output to input, a fune- 


tion of angular frequency, w 
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Fic. 2. Arbitrary transient input 
(PA)o; = phase angle between output and input, a 
function of angular frequency, w 
wo = 2x 7 = fundamental angular frequency 
— PA . ° ‘ 
(TF)o, = (AR)o, e’*’! = transfer function, a 


function of angular frequency, w 


When the input to a device or system is given, it can 


be represented as in Eq. (1) by 


i= DY Jae (3) 
n = 
where 
a 
Ia) = (1/T) / Ie" di (4) 
for 2s = 0, #1, #2, +Z, «....... 
Since the transfer function may be written 


(AR) je? *, the output component for the nw» fre- 


quency will be 


J( n yet nol 4 R)o(n el! PA yg) 
J(n)(AR)o(nye 


i{nant+(PA),7(")] 


which accomplishes the required amplitude multiplica- 
The frequency w in question is for 
therefore, the transfer function is 


tion and phase shift. 
the nth harmonic; 

identified by the coefficient ” as well as w (w = 
Now, summing the terms for all input frequencies, as 
in Eq. (3), the total output is 


Nw). 


JI(n)(AR)o (net FA! (5) 
where the values of /(m) are determined from the in- 
put as per Eq. (4). 

The theory for the inverse situation (second case) 
is as follows: Let the output from a dynamic system 
and the input to it be known, and let it be required that 
the overall frequency response (i.e., the transfer func- 
tion) be determined. 

The output may be written as a series just as the 
input was in Eqs. (3) and (4): 


n= + 


> P(e" (6) 


where 
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T 
P(n) = (1 nf O(t)e'"*" dt (7) 
0 


and 7 is the same period as in Eq. (4). But, since the 
output is also given by Eq. (5), we have 


n= +o , 
[J(n) (AR)o)(n)e\?4 ol) | inet 


] ol , 
P(n)e"™ (8) 
n= — o@ 
Now, since 7, the period, is the same in both Eq. (5) 
and Eq. (7), the fundamental frequency, wo, is, of course, 
the same in both. 


Since Eq. (8) must hold for all values of the variable, 
the two series must be equal term by term. ‘Thus, for 
any , 


J(n) (AR)o(nyb??o™ = P(n) 


and, referring to Eqs. (4) and (7), we have 
(TF)o(n) = (AR)or(n)e? 1 
= P(n)/J(n) 


» a 
= ¥ O(t)e mata] f I(t)he "dt = (9) 
0 / 0 


The transfer function at the frequency nw, is thus 
determined for any n. The (7°F)o,; for a frequency be- 
tween two integer multiples of the originally selected 
fundamental frequency, wo, can be calculated, if it is 
required, by using a slightly different wo. It should be 
emphasized that the period, 7, must be the same in 
both expansions, since Eq. (9) depends on this condi- 
tion. The relative accuracy of the calculated value of 
the transfer function for a particular frequency will, of 
course, depend on the relative contribution of the terms 
of this frequency in the known input and output. 


A convenient and practical tabulating machine 
method for solving these problems is given in reference 
12, which makes use of IBM-matrix multiplication 


techniques for the numerical integration of Eq. (9). 

The transient response thus obtained provides in- 
formation on degree of stability, response time, accur- 
acy, etc.; and, therefore, it can be seen that the Fourier 
technique plays an important role in the application of 
the TFF method to the analysis of a dynamic system. 

The TFF method in general, then, provides both the 
forced harmonic response of the system and the tran- 











8, MINUS 9, _ 80, 4 
B 
(INPUT) A (OUTPUT) 
Fic. 3. Block diagram of simple feedback servomechanism with 


frequency dependent elements in forward and feedback paths. 
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sient response to arbitrary inputs, which results in 
complete solution to the physical problem. 

A brief summary of the history leading to the devel- 
opment of the TFF method in the field of automatic 
control will now be given. 


(c) History 


The beginnings of closed loop analysis can be traced 
back to the basic work of Nyquist! and Black,’ who in 
1932 and 1933 applied these ideas to the problem of 
stabilizing feedback amplifiers. Their work eventually 
resulted in mathematical criteria for determining the 
dynamic stability of closed loop feedback systems in 
general. 

Much of the original work in the field of servomech- 
anisms, which was an outgrowth of automatic con- 
trol analysis, was done by Hazen* in 1934; in fact, it 
was Hazen who first introduced the term ‘‘Servomech- 
anism.’’ However, Taplin,‘ in 1939, was one of the 
first to notice the similarity between the analysis of 
servomechanisms and feedback amplifiers. 

With the beginning of World War II, the interest in 
servomechanisms and automatic control was greatly in- 
creased, with the result that highly skilled scientific per- 
sonnel were called upon to develop further the basic 
theories. Among those active in this further develop- 
ment were Bush, Hazen, Gardner, Draper, Brown, 
Harris, Hall, Campbell, Kochenberger, Taplin, Phil- 
brick, Adrendt, Nyquist, Bode, Black, Ferrell, Dietzold, 
MacColl, and Minorsky, to name a representative num- 
ber of the outstanding contributors to the art. A 
complete bibliography covering the work of these men 
is given in reference 5. 

Following World War II, the basic concepts of auto- 
matic control and dynamic response of multiple-loop 
systems were extended further, especially in their ap- 
plication to the automatic control of aircraft. Much 
of this work was done by W. F. Milliken, Jr., E. V. 
Laitone, R. C. Seamans, B. G. Bromberg, and others.®~* 

The Transfer Function-Fourier method, which was 
an outgrowth of this basic research, was then made 
practical by the further development of graphical and 
methods for handling the Fourier tech- 
12 


machine 
nique.” 

Before discussing the extension of the TFF method 
to the dynamic analysis of aeroelastic aircraft, a few 
examples of the conventional application of this method 
to the analysis of simple servomechanisms and auto- 
matic control systems will be given for comparison later 


in this paper. 


CONVENTIONAL APPLICATION OF TFF METHOD 
TO SIMPLE FEEDBACK SERVOMECHANISMS 
AND AUTOMATIC CONTROL SYSTEMS 


(IIT) 


The conventional application of the TFF method 
to problems in automatic control can be demonstrated 
by first considering the case of a simple feedback servo- 
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mechanism with frequency dependent elements in 
both the forward and feedback paths. Fig. 3 illustrates 
this case. 

Note: Subtraction is usually obtained at point A 
by changing the sign of 8 and adding—.e., the phase 
angle of 8 is shifted by 180°. 


6; = input 

A = output 

m = forward transfer function 

B = feedback transfer function 

up = open-loop transfer function 

u/(1 + uB) = 6/0; = closed-loop transfer function 


The overall (closed-loop) transfer function for the 
system may thus be computed directly from the ele- 
ment transfer functions, u and 8, using the conventional 
Nyquist methods* for checking the stability of the sys- 
tem. From then on,this system and its overall trans- 
fer function may be considered as an element or sub- 
system in a larger system. 

A typical example of a larger system for automatic 
control of a command guidance aircraft is shown in 
Fig. 4. 

In this case the input to the system is the commanded 
acceleration, ”,, and the output is the actual acceleration 
response, %, at the aircraftc.g. The individual transfer 
function for each element in the system may be obtained 
by the corresponding ‘‘expert’’ who knows most about 
that element. For example, (1) the transfer function 
of the amplifier plus actuator combination, including 
simulated control surface load feedback, may be ob- 
tained directly from laboratory tests without resorting 
to theory; (2) the aerodynamic transfer functions be- 
tween harmonic control surface motion and aircraft 
response (6, 6, mo, etc.) may be obtained from the con- 
ventional (see reference 9, Appendix H) differential 
equations for aircraft motion; (3) the instrumentation 
transfer functions for the accelerometer, the rate gyro, 
and the control surface pickoff may be computed di- 
rectly from damping ratio and natural undamped fre- 
quency specifications for these instruments (with the 
knowledge of the impressed voltages across each pick- 
off); and (4) the electrical shaping transfer functions 
may be computed directly from knowledge of the elec- 
trical networks and their constants. 

The overall (net closed-loop) transfer function of the 
system may then be computed from the individual 
transfer functions, using the usual Nyquist techniques 
for a multiple-loop servomechanism,t which include 


* Nyquist showed that when the open-loop (uf) transfer func- 
tion for this simple system is equal to —1, or encircles the point 
[—1, 7(0)] in the Nyquist plot, the closed-loop transfer function 
becomes unstable. 2 17 

+t The Nyquist criterion for stability of a multiple-loop direct 
feedback servomechanis:n is that the net number of encirclements 
of the point [—1, 7(0)], considering progressively the contribution 


of each open-loop transfer function, be equal to zero." * 17 
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Fic. 4. Block diagram for command guidance of aircraft. 
m = command acceleration; m) = output acceleration at c.g.; 
Na = acceleration at point a, / ft. ahead of c.g.; E. = command 


shaped command voltage; Ep = control surface 


voltage; E., = 
pickoff voltage; Eps = shaped control surface pickoff voltage; 


Ee = gyro voltage; Egs = shaped gyro voltage; E4 = accelerom- 
eter voltage; Eas = shaped accelerometer voltage; E = volt- 
age from amplifier to actuator; 6 = control surface displacement ; 
6 = pitch angle of aircraft; 6 = rate of pitch; and 6 = accelera- 
tion in pitch. 


checking the stability of the “‘open-loop”’ transfer func- 
tions as each successive loop is added to the system. 

The overall transfer function of this larger system 
may then be used with the Fourier technique to com- 
pute the transient acceleration response of the aircraft 
to an arbitrary command or, it may be used as an ele- 
ment or subsystem in an even larger automatic remote 
control system, which might include the degrees of 
freedom of a radar system, a computer, etc. This 
ability to treat each part of the system separately and 
to consider the effect of successively adding each part 
is one of the primary advantages of the TFF method. 

In summary, these two cases—the simple servomech- 
anism and the larger system for automatic control of 
a command guidance aircraft—demonstrate the con- 
ventional application of the TFF method. 

It will now be shown that, by extending these funda- 
mental concepts, a rigid-body and an aeroelastic air- 
craft may be analyzed dynamically by the TFF method, 
and further, that an aircraft stability analysis, a dy- 
namic tail-load analysis, and a flutter analysis may be 
obtained in only one analysis. Also, from the foregoing 
discussion, it can easily be seen that it is especially con- 
venient in this application of the TFF method to treat 
the aircraft itself as a servomechanism. 
(IV) Extension oF TFF Metuop To DyNaAmIc 
ANALYSIS OF RIGID-BODY AND AEROELASTIC AIRCRAFT 

(TREATING THE AIRCRAFT ITSELF AS A 
SERVOMECHANISM) 


One does not usually think of the aircraft itself as a 
servomechanism. However, it is convenient to do 
this in applying the TFF method. 
on an airplane cause it to move, and the motion in turn 
“feeds back’”’ to affect the forces; or, forces on a wing 
cause wing deflection and airplane motion—but wing 
deflection ‘‘feeds back”’ to affect the forces on the wing 


and therefore the airplane motion, and the airplane 


For example, forces 
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Fic. 5. Block diagram for horizontal tail loads produced by 
elevator deflections (neglecting aerodynamic circulation lag and 
aeroelastic deflections). 
motion ‘‘feeds back’’ to affect the forces—and, there- 
fore, the block diagram for the system is basically no 
different from that for a feedback servomechanism or 
an automatic control system. Therefore, by extending 
the fundamental concepts used by the TFF method, a 
rigid-body or an aeroelastic aircraft may be analyzed 
dynamically by treating it as a multiple-loop servo- 
mechanism with aerodynamic and aeroelastic feedbacks. 

In this way, many of the methods and ‘“‘tools’’ that 
were developed years ago by the electrical-servomech- 
anism engineers in their study of the dynamic stability 
of communication networks, feedback amplifiers, servo- 
mechanisms, and automatic control systems may be 
applied directly by the TFF method to the dynamic 
analysis of rigid-body and aeroelastic aircraft. 

Also, as shown in Part (IV-c) of this paper, many of 
the analysis techniques recently developed in the study 
of the dynamics of aeroelastic aircraft by conventional 
methods, especially the matrix techniques developed 
in the field of flutter, may be applied in the dynamic 
analysis of systems for automatic control. 

A simple rigid-body application will now be given to 
demonstrate the extension of the TFF method: 

(a) Example of the application of the T FF method toa 
rigid-body longitudinal dynamic analysis of an airplane 
to determine the transient maneuvering tail loads devel- 
oped by elevator displacements. 

A simplified block diagram for the tail loads pro- 
duced by elevator deflections is shown in Fig. 5. 

6. = elevator displacement 


6 = rate of pitch of airplane 


Uy = initial velocity of airplane 

AW = normal velocity of airplane, w = AW/U% 

€ = downwash angle 

a, = tail angle of attack, referred to zero lift line 

be = tail length 

7 = horizontal tail load 

n = tail efficiency = dynamic pressure at tail + 
free-stream dynamic pressure 

Cz = normal force coefficient 


Note that elevator motion is the input to the system 
and horizontal tail load is the output. Also, in this 


simplified case the effects of aerodynamic circulation 
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lag and aeroelastic deflections have been neglected 
i.e., it has been assumed that the aircraft stability de- 
rivatives are independent of frequency as well as amphi- 
tude. This last assumption allows the use of the con- 
ventional (see reference 9, Appendix H) aircraft differ- 
ential equations in computing the airplane transfer 
functions between elevator deflection input and airplane 
Also note 
that this simplified block diagram treats the system as 
an “‘open-loop”’ 


normal velocity and rate of pitch outputs. 


i.e., without feed- 
backs (actual feedbacks are included in the differential 
equations for the transfer functions of the airplane as a 
whole). 


servomechanism 


The basic equations relating the transfer functions 
for the system are: 


— feCs\ 1 . eo 
(TF)ars, = ( Ow ) 2 pS, Uo'n( TF) seus, 
where 
CF il . ‘T° ~ l 
(TF) naa = (TF ws, _ (*) x 
da 


sittin l, ) an | l, oe") 
Pia, ~ (= krmn,|+ —- (rp, ( 
| . (; ere tare ne Tos. 


Then, once (TF), ,, (TF) 
puted from the conventional aircraft differential equa- 


» and (7F)js, are com- 
tions, it is a simple matter to compute (7°F)47s,, which 
is then all that is needed to determine the transient AT 


° 


to any 6, (using the Fourier technique). In this way, 
the differential equations need only be solved once to 
obtain the transfer functions—i.e., every time the in- 
put, 6,, is changed they do not have to be solved over 
again to obtain the corresponding AT. 


The TFF method for computing horizontal maneu- 
vering tail loads has recently been applied by the Doug- 
las Aircraft Company to practical examples, and the 
method has been compared with the older method out- 
lined in A.A.F. Technical Report No. 5185. The com- 
parison showed that the results of the two methods 
agreed to within 5 per cent. The comparative ad- 
vantages, however, of the TFF method are the time 
savings and the ability to handle an arbitrary input. 
The A.A.F. method is only conveniently applicable to 
trapezoidal type elevator inputs, and the system equa- 
tions must be resolved each time the input is changed. 
(IBM calculating time for the TFF method usually 
amounts to about one hour per input. On the other 
hand, the time required per input for the A.A.F 
method is approximately 40 hours for the first input, 
with approximately 5 hours for each additional input.) 

This, then, demonstrates only one of the many appli- 
cations of the TFF method to a rigid-body analysis. 

A rigid-body analysis may easily be extended, how- 
ever, to include the degrees of freedom of aerodynamic 
circulation and aeroelastic deflections; a simple ex- 
ample of this extension of the TFF method is as follows 








th 
est 
wi 
sp: 
fac 
be 
th 
for 
ize 


flu 


ins 
be 


thi 
tic 
as 
pu 
da 


fli; 
me 
th 
gi) 
in 


no 
in 


“ted- 
ty de- 


umpli- 
eC con- 
differ- 
ansfer 
‘plane 
» note 
em as 
feed- 
ential 
Pasa 


‘tions 


4 
06, 
‘om- 
qua- 
hich 
AT 
vay, 
e to 
in- 
ver 














DYNAMIC ANALYSIS OF 


(b) Example of the application of the TFF method to 
a longitudinal dynamic stability and performance anal- 
ysis (and a limited flutter analysis) of an airplane 
equipped with a mechanical gust load alleviator that is 
actuated by elastic wing deflection. 

A simplified block diagram for the entire gust load 
alleviation system is shown in Fig. 6, where the 
rigid-body feedbacks, 2,, my, 2, and m,, have been 
included in the element for the airplane as a whole. 
The corresponding system transfer function equations 
are given on page 388. In this case, for simplicity, the 
transfer functions are denoted by ‘“‘F”’ instead of “TF.” 
The operation of this simple device is as follows: The 
lift caused by the gust creates an elastic deflection of the 
wing which in turn moves the ailerons to reduce the lift 
In this sense, the alleviator itself 
Since the opera- 


caused by the gust. 

acts as a ‘‘feedback”’ in the system. 
tion of the device depends upon wing flexibilitv, any 
analysis of its behavier must consider not only the mo- 
tion of the aircraft as a whole, but the important aero- 
elastic effects as well. Also, since the air speeds are 
comparatively low and the aeroelastic frequencies are 
comparatively high, the effects of the aerodynamic 
(circulation lag) transfer functions must be considered. 
Note that the block diagram for this system is of the 
closed multiple-loop feedback type and that it is basi- 
cally no different from that for the automatic control 
system discussed in Part (III) of this paper (see Fig. 4). 

The input to the system is either the gust velocity or 
the elevator motion, and the outputs of primary inter- 
est are the acceleration, , at the airplane c.g. and the 
wing (first mode) bending deflection, d, at a particular 
spanwise location. In considering wing deflection, the 
facts that there are spanwise distributions and that both 
bending of, and torsion about, the elastic axis enter into 
the first bending mode, make it convenient to account 
for this by using a generalized force, P,,, and a general- 
ized wing deflection, d, which are commonly used in 
flutter analyses. '*—™® 

The analysis of this system by the TFF method, treat- 
ing it as a closed multiple-loop servomechanism, has 
been completed for the C-47 airplane by the Douglas 
Aircraft Company, and a flight test program to test 
the system is in progress. Elevator motion, in addi- 
tion to, but separate from, gust velocity, was considered 
as an input in the analysis to provide theoretically com- 
puted calibration and simulated evaluation response 
data to be used for comparison later with corresponding 
flight-test response data for elevator inputs. 

An explanation of the transfer functions for each ele- 
ment in the block diagram, with information concerning 
the method of computing each transfer function, is 
given in Appendix A of this paper. Also included in 
Appendix A is a table of symbols and conventions used 
in the gust-load alleviation analysis. 

The actual transfer functions for each element are 
They are, however, included 
Each transfer function 


not given in this paper. 
in Appendix II of reference 16. 
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was combined vectorially in accordance with the sys- 
tem transfer function equations shown on page 388. 
This was done in the usual graphical manner to first 
obtain the open-loop and then the overall (closed-loop) 
transfer functions. 


Open-Loop Transfer Functions for Gust Inputs 


Figs. 7 and 8 show the open-loop transfer functions 
for Fig and Fi» at indicated air speeds of 150 and 190 
m.p.h., respectively. Referring to the system transfer 
function equations on page 388, the open-loop transfer 
function for both Fio and Fx) would be [FoF 34( Fs, + 
FesF 33) |, where the primary part is Fy FyFs. As the 
equations are written, the point [+1, 7(0)| would be 
the danger point on the Nyquist plot. However, to be 
consistent with the usual servomechanism convention, 
which is based on the point [—1, 7(0)], the negative of 
this open-loop transfer function was plotted in Figs. 
7 and 8. The open-loop plots for Fy and for d/ a, and 
n/a, at 150 and 190 m.p.h. (indicated air speed) showed 
no tendency at all toward encirclement of —1, 0 and, 
therefore, are not presented. None of the open-loop 
plots for an indicated air speed of 110 m.p.h. showed 
any tendency toward instability; therefore, they are 
not presented either. In all the cases shown on Figs. 7 
and 8, except the one specifically indicated on Fig. 8, 
the open loops were stable, as indicated by the usual 
Nyquist criterion,' of nonencirclement of the point 
— 1,0. 

Particular reference should be made to Fig. 8 showing 
the stability plot (open loop) for the Fy) and F2» transfer 
functions at 190 m.p.h. (IAS). A more complete 
investigation of the effect of change in the gear ratio 
between wing deflection and aileron movement was 
made at this speed since it was considered to be the most 
critical. While one cannot base a decision of adequate 
damping entirely on how close the different curves 
approach —1, 0, it is evident that, as the gear ratio is 
increased, there is a definite tendency toward encircle- 
ment of —1, 0—i.e., instability. This trend, together 
with consideration of the corresponding overall transfer 
function and the transient behavior, establishes upper 
limits for the gear ratio. 

The alleviator design considered in the analysis has 


the following characteristics: 
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WING DEFLECTION,¢ 


The various gear ratios used in the analysis were as 
sumed to be linear and are identified as shown. Table 
| shows the actual values of the gear ratios and the me- 


chanical stiffness associated with each. As long as the 
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transfer functions for the Alleviator 


Mechanism, Appendix A. 


Simplified block diagram for a gust load alleviation system. 
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transient aileron deflection remains within the limit 
shown, the assumption of linearity is valid. This re- 
quires that the aileron response to typical gusts be in- 
vestigated as part of the analysis. (This was done later, 
and it was determined that the aileron response was 
within the linear limit for the gust magnitudes investi 


gated. ) 


Overall Transfer Functions for Gust Inputs 


Figs. 9, 10, and 11 plus 12 show the overall transfer 
functions relating wing deflection (stress) to gust inputs 
for indicated air speeds of 110, 150, and 190 m.p.h., re- 
The various alleviator gear ratios used are 
The lowering of the amplitude 


spectively. 
indicated on the curves. 
ratio curves in the low-frequency region as compared 
to the “‘alleviator out’’ curves is an indication of the 
stress alleviation to be expected in the transient response 
to an actual gust. However, while an increase in the 
gear ratio reduces the low-frequency portions of the 
curve, note that the peaks near 5 cycles per sec. are in- 
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TABLE | 
Value of g* Stiffness, A * 

é, Rad. Ft.Lbs. 

Gear Ratio e Pe. Rad. 6. 
Average 0.22 3,060 
Present maximum 0.30 2,240 
1.5 X present maximum 0.45 2,240 
2.0 X present maximum 0.60 2,240 
4.0 X present maximum 1.2 2,240 


* These are the values of g and K used in the equations for Fy 
and Fy, (see Appendix A of this paper ) 


creased. This indicates an increasingly sharp resonance 
of the wing-aileron loop near the wing’s natural fre- 
quency, which means that the loop is becoming less 
damped—i.e., approaching instability. This is di- 
rectly related to the increasingly close approach to the 
point — 1, 0 on the stability (open-loop) plot. 

Figs. 13, 14, and 15 show the overall transfer func- 
tions relating aircraft c.g. acceleration to gust inputs 
for indicated air speeds of 110, 150, and 190 m.p.h., 
respectively. These show the same lowering of the 
amplitude ratio curves in the low-frequency region, but 
the progressive peaking of the curves at the higher fre- 
quencies is not as pronounced. This is to be expected 
since the aircraft response to wing deflection at the 
higher frequencies would be small. 


Transient Responses to Gusts 


As previously noted, the lowering of the overall trans- 
fer function curves in the low-frequency region, as 
compared to the case with the alleviator out, is an indi- 
cation of the alleviation to be expected in the transient 
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response to an actual gust. Furthermore, the serious- 
ness of the peaks in the wing deflection curves at the 


Clearly 


what is needed is the transient response of the system, 


higher frequencies is not immediately evident. 


both in wing deflection and airplane acceleration, to 
typical gust inputs. Comparison of these transients 
for various alleviator gear ratios with that for the case 
when the alleviator is disengaged will yield the per cent 
alleviation and the duration and magnitude of any os- 
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cillations at the higher frequencies. Consideration of 
the interaction of these effects will establish the most 


desirable gear ratio consistent with the compromise be 


bility. 

The method used to obtain the transient response 
was the Fourier technique discussed in Part (II-b) of 
this paper. In the application of this technique, use 
was made of the IBM-matrix method for numerical 
integration discussed in reference 12. 

Figs. 16 and 17 show the wing deflection transients in 
response to a sharp-edged gust, as shown in the follow- 








ing sketch, at air speeds of 110 and 190 m.p.h., respec- 


tively. 
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TIME ——o 





Sharp-edged gust. 


It has recently been discovered that the zero fre- 
quency values of the transfer functions F), and F 2; are in 
error. As a result, the overall transfer functions for 
d/ay and n/ay do not have an amplitude ratio of zero 
at zero frequency, as they would if correct values had 
been used. It has been proved that the error does not 
appreciably affect the per cent alleviation indicated on 
transient responses of Figs. 16-21 but causes the 
transients to fail to return to the initial steady state. 

In all of the figures, reference to the case of “‘alle- 
viator, out”’ refers to an airplane without an alleviator 
present. However, in the current C-47 flight tests, 
the alleviator will be present even though inoper- 
ative. Its mere presence changes the stiffness of the 
aileron system and changes slightly the dynamics of the 
wing response. In order to evaluate this change, the 
wing response was evaluated for the case of the allevi- 
ator present but inoperative. Fig. 17 shows that at 
190 m.p.h. the expected alleviation over an airplane 
with no alleviator is 15 per cent, while in comparison 
with an airplane having an alleviator present but in- 
operative, the alleviation is 19 per cent. This indi- 
cates that data to be obtained on the current C-47 
flight tests will show a little more alleviation than is 
actually present. At 110 m.p.h., as shown in Fig. 16, 
the alleviation was only 10 per cent as compared to 
15 per cent at 190 m.p.h. Both Figs. 16 and 17 show 
the tendency of the wing oscillations to last longer as 
the alleviator gear ratio is increased. 

The c.g. acceleration for the same sharp-edged gust 
is shown in Figs. 18 and 19 for air speeds of 110 and 
190 m.p.h., 
mation yielded by these curves is that the c.g. acceler- 
ation alleviation is less than the deflection (stress) allevi- 


respectively. The most important infor- 


ation at any air speed—i.e., 9 per cent for acceleration 
and 15 per cent for deflection at 190 m.p.h. 
cated in Fig. 18, the computed acceleration alleviation 
at 110 m.p.h. was too small to be considered. 

Wing deflection transients to a 30 ft. per sec. gust with 
the gradient shown in the following sketch were com- 
puted and are shown in Fig. 
speed of 190 m.p.h. 


As indi- 


20 for an indicated air 
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Transients were obtained in this case for a large num- 
ber of different gear ratios to help establish the maxi- 
mum allowable gear ratio consistent with adequate 
system stability. As shown in Fig. 20, alleviations of 
16, 20, 27, 32, and 37 per cent were obtained for the 
average gear ratio; present maximum gear ratio; and 
one and one-half, two, and four times the present maxi- 
mum gear ratio, respectively. However, accompanying 
this increasing alleviation is an increase in the oscillation 
time of the wing. (This increased oscillation time is 
directly associated with the peaking of the amplitude 
ratio curves in the overall transfer function for wing 
deflection, Fig. 11, near 5 cycles per sec.) In fact, a 
gear ratio of four times the present maximum was de- 
liberately chosen to show an unstable case. This case 
yielded an alleviation of 37 per cent, but note that the 
oscillations did not die out but slowly built up. This 
would mean wing flutter due to coupling of the allevi- 
ator with the first bending mode. Obviously, the gear 
ratio of four times the present maximum is too large. 
In fact, the gear ratio of twice the present maximum 
causes oscillations that last an undesirably long time. 
Those for 150 per cent of the present maximum do not 
last too much longer than those associated with the 
present maximum. Reference again to the amplitude 
ratio curves in the overall transfer function for wing de- 
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flection, Fig. 11, and the shape of the peak associated 
with the 150 per cent curve, indicates that it is in the 
gear ratio region where a rapid change in peak ampli- 
tude may be expected with further increase in gear 
ratio. Therefore, it would seem reasonable to con- 
clude that the gear ratio near zero deflection could be 
increased to one and one-half times the present maxi- 
mum ratio and still maintain a reasonable stability 
margin. The gain in deflection (stress) alleviation 
from 20 to 27 per cent appears to warrant such a change. 
However, the slow gain in alleviation from there on, 
coupled with the danger of making the system un- 
stable, does not warrant further increase in the gear 
ratio. 

The c.g. acceleration transients in response to the 
same 30 ft. per sec. gradient gust at 190 m.p.h. are shown 
in Fig. 21. Only the cases for the present average gear 
ratio and the 50 per cent increase in present maximum 
gear ratio are shown. As previously found for the case 
of a sharp-edged gust, the acceleration alleviation was 
smaller than the deflection alleviation (15 per cent ac- 
celeration alleviation as compared to 27 per cent deflec- 
tion alleviation for the case of the 50 per cent greater 
gear ratio). 

In summary, the analysis showed that the system to 
be tested was stable and that it was more of a wing stress 
alleviator than an aircraft acceleration alleviator. It 
also showed that the alleviator gain could be increased 
to one and one-half times the present maximum gear 
ratio and still maintain a reasonable stability margin, 
and further, that if the gear ratio were increased by 
two and one-half times the present maximum value, 
the alleviator plus wing combination would flutter (at 
about 5 cycles per sec.) and any further increase would 
In addition, 
the allevi- 

The suc- 


cause the system to become unstable. 
the analysis gave the performance of 

ator (i.e., the per cent alleviation, on-to-off). 
cess of this analysis also showed that the TFF method 
is a practical way of handling dynamic aeroelastic prob- 


lems. 


Equations for Gust-Load Alleviation System 


F, = FF is + F Fas Fx = FyFy4 + F Fog 
Fe = Fik\g + Fikes Fy = FyFi3 + F5 Fo: 
Fe = 14+ (FoF + Fi Fx) Fy = F3Fs + FiFe 
Fp = FyFy + PsP Fy = F\Fy + FoFs 
Fy = Fax/(1 — Fa7F33) 
— ba vs FoF + FoPs( Ps + Fos Fis) 
a hin [FoF ss( Fa + Fos Fss) | 
—— d > Fo( Fi + Fu Fas Fos) 
e w+ ay 1 — [ Fo Fsa( Fs + Fos Fs3) | 
for 6, = (), 
d Foo Fe 


ay 1 — (Foo F'4 + FioF z) 
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n = Fp 4 rites + Fils) 
ay 1 — (Foo F's + Fy Fp) 
for ay = 0, 

d FoF n 

be ie (Foo F 4 + Fy Fz) 

n Fy( Foo Fy FoF xz) 
‘eas F, r H( = as = = 

be l= (Foo F 4 + FoF gz) 

Although the analysis of the gust-load alleviation 


system showed that the alleviator plus wing combina- 
tion could be made to flutter, it was felt that this was 
by no means a complete flutter analysis because higher 
frequency elastic modes had been neglected. How- 
ever, this led to the idea that it would be possible to 
include these additional elastic modes in future anal- 
yses by this method. This then prompted the idea 
of investigating the possibility of using the TFF method 
to obtain, in a single analysis, an aircraft stability 
analysis, a dynamic tail-load analysis, and a_ flutter 
analysis. 

For the past year and one-half this possibility has been 
investigated at the Douglas Aircraft Company, and it 
has been found to be entirely feasible. An example 
of this further extension of the TFF method is as fol- 
lows. 

(c) Example of the application of the T FF method 
to a longitudinal dynamic analysis of an aircraft, to de- 
termine by a single analysis (1) the stability of the air- 
craft flight path, (2) the transient tatl loads, and (3) the 


flutter response of the control system and the multielastu 


structure, where all three of these responses result from an 
arbitrary gust disturbance or a control force input. 

The investigation of the possibility of using the TFF 
method to make this more comprehensive analysis was 
stimulated further by the need for this analysis. The 
need can be explained as follows: A rigid-body free- 
control longitudinal dynamic stability analysis for one 
of the Douglas airplanes had previously been made, 
including the degrees of freedom of the airplane, the 
elevator, and the elevator control tab systems. This 
analysis assumed that the aircraft structure was rigid 
and neglected the effect of aerodynamic circulation lag. 
It was found, however, that some of the natural frequen- 
cies obtained from this rigid-body analysis coincided 
with the aeroelastic frequencies of the structure. This, 
of course, made the results of the rigid-body analysis 
questionable. Simultaneously, it was discovered by 
Smilg* and others that conventional high-frequency 
flutter analyses could be affected appreciably by in- 
cluding the aircraft body modes. These two basic 
reasons therefore established the need for a more com- 
prehensive analysis that would extend the rigid-body 
analysis to include the high-frequency aeroelastic and 
circulation lag effects or, looking at it from the other 


* Chief of the Dynamics Branch, Aircraft Laboratory, A.M.C., 
Wright-Patterson Air Force Base, Dayton, Ohio 
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point of view, would extend the flutter analysis to in- 
clude the effects of the low-frequency body modes. 

It was then suggested by J. Kelley, Jr.* that the 
comprehensive analysis might also determine the ef- 
fect of elastic deflections and circulation lag on transient 
tail loads. 

This ‘overlap’? problem between high- and low-fre- 
quency response is, of course, much more pronounced 
in a large aircraft in which a greater proportion of its 
mass is distributed out along its wings. This results 
in lower aeroelastic frequencies and larger elastic de- 
flections. This trend is also being observed in the de- 
sign of flexible wing jet bombers with wing-tip engine 
pods. 

A block diagram of the more complete system for 
the longitudinal case is shown in Fig. 22. This dia- 
gram represents the airplane as a large servomechanism 
with aeroelastic, aerodynamic, and inertia coupling 
feedbacks, and, as can easily be seen, it is basically the 
same type of diagram as that for the automatic control 
system shown in Fig. 4. It is also similar to that for a 
feedback amplifier. 

The individual transfer functions for each element 
in the block diagram of Fig. 22 are explained in Appen- 
dix B, in which is also included a list of the symbols and 
conventions used in the analysis. It will be noted 
that many of the obviously unimportant feedbacks are 
not included in the block diagram, as for example Fy, 
Fa, Fs, Fos, Fos, etc. However, these could easily be 
included if required. The degrees of freedom accounted 
for are those for the aircraft itself, its elastic wing (first 
bending and first torsional modes), fuselage flexing 
(first bending mode), horizontal tail deflection (first 
bending mode), symmetrical aileron deflection, ele- 
vator deflection, elevator control tab deflection, aero- 
dynamic circulation lag, and control system dynamics. 
Also included are the effects of downwash lag, inertia 
coupling, and wing-to-tail lag in entering a gust. 

The basic inputs to the system shown in Fig. 22 are 
gust velocity or stick force, it being assumed that in a 
practical case only one input need be considered at a 
time. However, both inputs could be considered simul- 
taneously if required. 

The outputs from the system fall into three basic 
categories. These are as follows: (1) outputs asso- 
ciated with the low-frequency response of the aircraft 
as a whole—i.e., those which describe the stability of 
the aircraft flight path; (2) outputs associated with the 
transient response of the tail; and (3) outputs associ- 
ated with the flutter response of the control system and 
the multielastic structure. 

In considering the flight-path stability, the outputs 
(in the first category) would be as follows: The low- 
frequency flight-path response the 


1.€., Phugoid 


* Chief of the Design Criteria Unit of the Structures Branch, 
Aircraft Laboratory, A.M.C., Wright-Patterson Air Force Base, 
Dayton, Ohio. 
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would be dominant in the output of change in speed: 
while the higher frequency response—i.e., the 
would be dominant in the output of air- 


“short 
period’ mode 
craft c.g. normal acceleration. 
the second category, the transient tail loads, tail acceler- 
ations, and tail deflections (including fuselage and hori- 
zontal tail bending) would be outputs valuable in struc- 
tural design. The third category—i.e., the high-fre- 
quency flutter response of the system 
sented by outputs such as wing, tail, and fuselage de- 
flections, and aileron, elevator, and elevator control 


Likewise, in considering 


would be repre- 


tab deflections, etc. 

A practical example of an input for which the out- 
puts in these three categories are of interest is the con- 
ventional “‘short period’’ stability check used in proto- 
type flight testing. In this check, the stick force input 
is applied quickly to move the elevator from the trimmed 
flight position to a new position in about 0.2 sec., fol- 
lowed by sudden release of the stick. The resulting air- 
craft stability, the transient tail loads, and the flutter 
response of the system, are then, of course, the primary 
outputs of practical interest. 

Another practical example would be the case of an 
airplane suddenly penetrating a sharp-edged gust. 
The resulting flight-path stability, the transient tail 
loads, and the flutter response, would be of prime in- 
terest. 

Now, since the outputs from these three categories 
occur simultaneously, it would seem logical to compute 
the outputs simultaneously and thus avoid the ‘‘over- 
lap” problem. 

The basic block diagram shown in Fig. 22 can be 
either expanded for further detail or reduced. For ex- 
ample, the expansions of Fy; and Fy, are shown in the 
lower left-hand corner of Fig. 22. Likewise, the ex- 
pansion of F;; (not shown) would show the effect of 
aerodynamic circulation lag as well as inertia coupling. 
Or, the expansion of F2; (also not shown) would show 
the contributions to \/, of the wing, the fuselage, and 
the tail (including the effect of wing-to-tail time lag 
in entering a gust). On the other hand, the basic block 
diagram of Fig. 22 can be reduced to a simpler block 
diagram (see Fig. 23) by consolidating the direct feed 
back loops. Also included in Fig. 23 are the equations 
that relate the transfer functions for the reduced sys- 
tem. 

Consistent with this first reduction, Nyquist plots 
can be made to check the stability of each of the direct 
feedback open-loop transfer functions (Ai, As, As, - . - )- 
Then, the corresponding closed-loop transfer functions 
(u41, M2, M3, . . . ) can be computed for use in the reduced 
diagram of Fig. 23. 

A convenient method for multiplying and dividing 
transfer functions is the Log Modulus Contour method 
given by Brown and Campbell (reference 5, page 256). 
By this method, multiplication and division in linear 
coordinates become addition and subtraction, respec- 


tively, in logarithmic coordinates. 
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Although the reduced block diagram of Fig. 23 is still 
rather complicated, it is, nevertheless, simpler than the 
first and can be made even simpler by consolidating the 
cross feedback loops. Figs. 24 and 25 show the cor- 
responding block diagrams that result from cross feed- 
back consolidation, for control force and gust inputs, 
respectively. A convenient IBM-matrix technique for 
computing the cross-feed open-loop transfer function 
can be applied. This technique is one that makes use 
of matrix and machine methods to solve the transfer 
function equations in Fig. 23 by the process indicated 
in Figs. 24 and 25. In this way, the cross-feed open- 
loop transfer function, \ = 1 — A, can be computed 
over the frequency range of interest (for use in the Ny- 
quist plot), and the final closed-loop transfer functions 
can then be obtained as indicated by the equations in 
Figs. 24 and 25. 

Once the overall transfer-function for each output 
(nine in this case) has been obtained, the corresponding 
transient response to an arbitrary input can be com- 
puted by the Fourier technique discussed in Part (II-b) 
of this paper. 

For the overall system to be stable, it is necessary 
that the final cross-feed open-loop transfer function, X, 
neither pass through nor encircle the point [+ 1, 7(0)] in 
the Nyquist plot. As mentioned previously, the point 
| —1, 2(0) | would be of interest if —\ were plotted. 
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It is possible for some of the direct feedback open- 
loop transfer functions (Ai, Ae, As, . . . ) to indicate in- 
stability, and yet the overall system could be stable. 
This condition would, of course, indicate a system that 
is not “‘fail safe.”’ 

It is also possible for some, or all, of the direct feed- 
back open-loop transfer functions (Aj, Ao, As, . . - ) to indi- 
cate stability, and yet the overall system could be un- 
stable as indicated by the final crossfeed open-loop 
transfer function, X. 

Therefore, by making a systematic check on the sta- 
bility of all of the open-loop transfer functions, any 
which are found to indicate instability may be picked 
from the system and examined to determine the cor- 
rection required for stability. 

In summary, then, the entire system is progressively 
simplified, while the open-loop Nyquist plots are made 
to systematically check the stability of the system. 

An example of the application of the TFF method to 
this more comprehensive analysis has recently been 
-arried out by the Douglas Aircraft Company with 
excellent results. 


(V) CONCLUSIONS AND RECOMMENDATIONS 


In conclusion, the author would like to summarize 
some of the advantages of the TFF method, especially 
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DYNAMIC ANALYSIS OF 


in its application to the dynamic analysis of rigid-body 
and aeroelastic aircraft: 

(1) It provides, in a single analysis, information 
that is usually obtained in many separate analyses. 
In this way it also eliminates the ‘‘overlap’’ problem 
between high- and low-frequency response—i.e., it in- 
cludes the effect of body-modes on flutter and vice 
versa. 

(2) It provides, through the block diagram and 
the transfer function coricept, a clearer insight to the 
true physical picture of the dynamic system. The 
conventional flutter analysis tends to be somewhat 
mathematical making the physical interpretation 
difficult. It is believed that the TFF method places 
flutter on a simple engineering basis, which allows the 
system behavior to be more easily visualized. 

(3) It makes use of well-established tools and 
techniques developed in the electrical and servomecha- 
nism fields years ago. * 

(4) It is readily adaptable to automatic computa- 
tion and system simulation; the block diagram itself, 
or its simpler derivatives, may be used directly in 
setting up a simulator. 

(5) It provides response data 
checked by dynamic flight-test techniques. 

(6) It provides response data in a form that is 
normally required for the analysis and synthesis of 
i.e., the aero- 


that are easily 


systems for automatic aircraft control 
elastic aircraft itself may be just one element in a larger 
system for automatic control, and in this way the aero- 
elastic effects on the automatic control problem are 
accounted for. 

(7) The method does not require ‘high powered”’ 
mathematicians to do the work. 

(8S) Individual transfer function computations can 
be done by people who are experts in their own special 
field, and work on each part of the system can proceed 
separately and simultaneously. 

(9) Any particular part of the system may be 
simulated while other parts, which have better experi- 
ment-to-theory justification, may be handled by theory 
alone. In many cases only laboratory transfer func- 
tions are available for a particular element in the sys- 
tem—i.e., there may be no simple analytical expression 
available for the transfer function for this element, in 
which case the graphical plots of the laboratory trans- 
fer functions may be used directly. 

(10) The method can handle conveniently an arbi- 
trary input, as long as this input can be represented by 
a Fourier equivalent. Also, an additional arbitrary 
input may be handled with little extension of time 
(about '/, hour on the IBM or a Fourier simulator), 
without having to resolve the system equations. 

(ll) The zero-frequency (static sensitivity) value 
of each transfer function that has been computed from 
dynamic theory can be readily checked with static test 
data; this is especially convenient when checking zero- 
frequency values of aerodynamic transfer functions 
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computed from flutter theory, with the values obtained 
from static wind-tunnel data. 

(12) The systematic and successive ‘‘open-loop”’ 
(u8) check on the stability of component loops gives a 
good insight to the “‘fail-safe’’ characteristics of the 
system, as well as a good background for ‘‘trouble- 
shooting’ the system in case the net transient response 
is unsatisfactory. 

(13) The actual time history of the net transient 
response of each output is obtained, thus providing 
quantitative data on (a) the degree of stability, (b) 
the response time, (c) the speed of response, (d) the 
accuracy, etc., which is far more than just the roots of 
the system equations as usually provided by matrix 
techniques alone. 

(14) Both matrix and graphical servomechanism 
techniques are used, rather than just one or the other. 
For example, the matrix techniques can be used to de- 
termine those individual transfer functions that are 
themseives complicated; or the matrix techniques may 
be used to determine the cross-feed open-loop transfer 
functions for the “‘us’’ Nyquist plot. In other words, 
the IBM-matrix methods can be used to reduce com- 
putation time when required, without loosing the ad- 
vantages of the open-loop Nyquist stability checks and 
physical interpretation provided by the servomecha- 
nism approach. 

(15) The TFF method gives both the forced fre- 
quency response of the system and the transient re- 
sponse to arbitrary inputs, which results in a complete 
solution to the physical problem. 

(16) Finally, although the time required to do the 
more complete aeroelastic analysis may be longer than 
the time required to do each individual analysis by con- 
ventional methods, the total time required for the larger 
analysis is probably less than that required for all of 
the individual analyses. Even if this were not the 
case, any additional time would certainly justify the 
advantages listed above. 

(It should also be pointed out that the TFF method 
may also be applied to the individual simpler analyses 
when the more complete analysis is not required. ) 

The author would also like to make the following 
recommendations: 

Further research and development of the TFF method 
in its application to rigid body and aeroelastic problems 
is reeommended. For example, basic work should be 
done on the lateral case as well as additional develop- 
ment of the longitudinal case. Also, other inputs 
should be considered, such as gun-blast loads, bomb 
release loads, landing loads, etc., and possibly additional 
outputs and additional degrees of freedom should be 
considered. 

With regard to the potentialities of developing fur- 
ther practical applications of the TFF method, the 
“the surface has hardly been 


author feels that 


scratched.” 
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DETAILS OF GusT-LOAD ALLEVIATION 
ANALYSIS 


APPENDIX A 


An explanation of the transfer functions for each ele- 
ment in the block diagram, Fig. 6, with information 
concerning the method of computing each transfer 
function, is as follows: 


F, = transfer function relating normal velocity of 
the aircraft as output to pitching moment 
as input 

F, = transfer function relating normal velocity of the 
aircraft as output to normal force as input 

F;, = transfer function relating normal acceleration 
of the aircraft (change in g’s) as output to 
pitching moment as input 

F,; = transfer function relating normal acceleration 
of the aircraft (change in g’s) as output to 
normal force as input 
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These four transfer functions (Fi, Fo, F;, and F,) 
were computed from the conventional differential equa- 
tions for aircraft motion (see reference 9, Appendix 
H), using the stability derivatives for the C-47 air- 
plane. 


F,; = transfer function relating pitching moment on 
the aircraft as output to elevator displace- 
ment as input 

F, = transfer function relating normal force on the 
aircraft as output to elevator displacement 
as input 


F, and Fs were calculated by the method presented 
in reference 15, using the M.A.C. of the horizontal tail. 
They were then corrected to fit the corresponding zero- 
frequency (static sensitivity) values for the stability 
derivatives, ms, and 25,, which had previously been de- 
termined from static wind-tunnel tests. A rather large 
correction was required to account for the small aspect 
ratio of the tail, wing downwash, and fuselage inter- 
ference. 


F\, = transfer function relating normal force on the 
aircraft as output to angle of attack from 
the gust as input 

F\; = transfer function relating normal force on the 

aircraft as output to deflection of the wing, 
first bending mode, as input 

transfer function relating normal force on the 


> 
I 


aircraft as output to aileron displacement as 
input 


F\, was first computed for the wings alone, with air- 
plane normal displacement as the independent vari- 
able.'% 
normal velocity, then angle of attack, the independent 
The normal force on the horizontal tail was 


This was then transformed to make airplane 


variable. 
computed in a similar manner but was corrected to 
allow for downwash lag, aspect ratio, wing-to-tail lag 
in entering the gust, and fuselage interference. The 
net normal force on the aircraft was found by taking the 
sum of the components from the wing and the tail. 
This sum was then corrected to fit the corresponding 
zero frequency data obtained from static wind-tunnel 
data. 

F\, was also corrected to agree with zero-frequency 
data. 

Fy, equals the transfer function relating pitching on 
the aircraft as output to gust angle of attack as input. 

Similar corrections were made to Fs; as were made to 
F\,, with an additional correction for the shift in axes 
from 35 per cent of the local chord, for which reference 
axis the original calculations were made, to the 22 per 
cent wing M.A.C. c.g. location, for which the overall 
analysis was carried out. Since there were no theo- 
retical calculations available for the frequency response 
of the moment contributed by the fuselage, the sum 
of the wing and fuselage moments was assumed to vary 
with frequency in the same manner as for the wing 
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DYNAMIC ANALYSIS OF 


The contribution of the wing, fuselage, and 


alone. 
tail was obtained for zero frequency and assumed to 
maintain the same ratio of magnitude throughout the 


frequency range. 

When the aircraft enters a gust, the tail is affected 
at a time after the wing. This effect was included by 
converting this time lag into a frequency dependent 
phase lag, which was then used as a correction to the 
tail transfer function before it was added to the transfer 
functions for the wing and the fuselage. 

Fx, equals the transfer function relating pitching 
moment on the aircraft as output to wing deflection, 
first bending mode, as input. 

This transfer function was also corrected for the shift 
of the reference axes from 35 per cent of the local chord 
to 22 per cent of the wing M.A.C. 

Fy, equals the transfer function relating pitching 
moment on the aircraft as output to aileron deflection 
as input. : 

This transfer function was also corrected to fit zero- 
frequency data on the corresponding stability derivative. 

F,, equals the transfer function relating generalized 
force on the wing, first bending mode, as output to 
angle of attack of the wing as input. 

As in the case for Fy, Fy, was first computed for air- 
plane normal displacement as an input, using references 
13 and 15, then changed to aircraft normal velocity, 
and, finally, wing angle of attack as input. 

F;,; equals the transfer function relating generalized 
force on the wing, first bending mode, as output to wing 
deflection, first bending mode, as input. 

This transfer function was calculated by the method 
of reference 14, taking into account the spanwise in- 


tegral effect. 


transfer function relating generalized force on 
the wing, first bending mode, as output to 


Fy = 


aileron displacement as input 
transfer function relating hinge moment on the 
aileron as output to wing angle of attack as 


Fy a 


input 
transfer function relating aileron hinge moment 
as output to wing deflection, first bending 
mode, as input 
transfer function relating aileron hinge mo- 
ment as output to aileron displacement as 
(Not shown in Fig. 6 because in- 
feedback, next 


Fy; _ 


Fy = 


input. 


cluded as an internal see 


column) 


Fy was also transformed to change the input from 
aircraft normal displacement to normal velocity, and, 
finally, to wing angle of attack; and Fy was also cor- 
rected to agree with the zero-frequency value of the 
aileron hinge moment stability derivative, Mahs,. 

Fy; equals the transfer function relating generalized 
wing displacement, first bending mode, as output to 
generalized force on the wing, first bending mode, as 


input. 


AEROELASTIC AIRCRAFT 


This transfer function was calculated from the cor- 
responding equation of motion of the wing with an 
inertia-damping-elastic degree of freedom in which 
the damping ratio was assumed to be 0.015 and the nat- 
ural undamped frequency /, = 4.67 cycles per sec. 
This damping ratio is equivalent to the value of the 
Theodorsen g = 0.03, usually assumed for structural 
damping in flutter analyses.'* 

F;, equals the transfer function relating aileron deflec- 
tion as output to wing deflection, first bending mode, as 
input. 

This transfer function accounts for the effect of the 
alleviator mechanism on aileron deflection and is essen- 
tially a “gear ratio,” although the mechanism dynamics 
were included. 

Fe, equals the transfer function relating aileron deflec- 
tion as output to aileron hinge moment as input. 

Fx, and Fe; were computed from the equation of mo- 
tion of the mass-elastic single degree of freedom system 
The 
aerodyiiamic hinge moment on the aileron resulting 
(Fy) is included here as an 


(for the alleviation mechanism) as shown below. 


deflection 
and is therefore not shown in the 


from aileron 
internal feedback 
block diagram of Fig. 6. 

Curves showing the actual] transfer functions for each 
element are included in Appendix II of reference 16. 


Equations for Alleviator Mechanism 
Total hinge moment on aileron = J/ = /,6,. 


M = 


4e2e0e08 


(gd — 6.) k + my, + mde — F446, 





\em 


™ pander 


my = Fa W’ + Fad = ma, 


}+s 


used in block diagram of Fig. 6. Thus, 


—I],w? 5, + Rig + Fu bg = gdk + m, — w*cdm, 


or 
gdk + Fyd + Fy W’ — w*cdm, 
a = ” i . 
- Iw : k T F's 
= dF, + ma Fes = 54, + 4a 
where 
F gk — cw*m, 5, 
"54 ans o “ = 
—I,w" + k + Fy, d 
l bn ba 
Fes = = 


—Tw +tk+ Fu Ma m,, (in Fig. 6) 


Basic Conditions Used in Analysis 


g = 0.220 rad. per ft. (average gear ratio) 


( = () 
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k 
Iq 


3,060 ft.lbs. per. rad. 
1 slug-ft.? 


Symbols and Conventions Used in Gust-Load Alleviation 
Analysis 


( = aileron unbalance, positive nose heavy, ft. 

d = generalized wing deflection, positive down (for first 
bending mode) 

Fy; = transfer function, output /input 

g = alleviator gear ratio (static sensitivity) 

he = effective moment of inertia of aileron (including 
alleviator mechanism) 

k = alleviator mechanism spring constarft 

m, = pitching moment on aircraft, positive airplane nose 
up 

n = aircraft normal c.g. acceleration, positive pull up 

P = generalized force on wing, positive down (for first 
bending mode) 

Uo = longitudinal component of aircraft velocity 

Vo = gust velocity, positive updraft 

W = aircraft normal velocity, positive down 

w = W/U> = nondimensional normal velocity, standard 
convention 

é ‘ = normal force on aircraft, positive down 

a = wing angle of attack, standard convention 

ba = aileron alleviation deflection, positive trailing edge 
down 

be = elevator deflection, positive trailing edge down 

Me = aileron mass unbalance, slugs 

Ww’ = W+ Vy 


APPENDIX B-—-DETAILS OF COMPREHENSIVE 
AEROELASTIC ANALYSIS 


Notation, Conventions, and Definitions 
The notation, conventions, and definitions of the transfer 
functions used in the comprehensive aeroelastic analysis are as 


follows: 


u = change in longitudinal velocity of airplane, 
positive forward 

0 = pitch angle of airplane, positive nose up 

w = W/Uo 

d, = elastic tail deflection (including first mode of 
longitudinal fuselage bending and first mode 
of horizontal tail bending), positive down 

d, = wing deflection for first bending mode, posi- 
tive down 

d, = wing deflection for second mode (first torsion), 
positive trailing edge down 

ba = aileron deflection, positive trailing edge down 

6, = elevator deflection, positive trailing edge down 

ber = elevator tab deflection, positive trailing edge 
down 

Xp = total longitudinal force on airplane, positive 
forward 

My = total pitching moment on airplane, positive 
nose up 

Z1 = total normal force on airplane, positive down 

i = net load on horizontal tail surface, positive 
down 

P, = generalized force on wing for first bending mode 

P» = generalized force on wing for second mode (first 
torsion), positive 2 causes positive dy 

Van = net moment at aileron hinge, positive aileron 
trailing edge down 

Ven = net moment at elevator hinge, positive ele 


vator trailing edge down 
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ll 


net moment at elevator control tab hinge, posi- 
tive elevator control tab trailing edge 
down 

aircraft normal acceleration at center of gray 
ity, positive upward (i.e., positive for pull 
up) 

aircraft normal acceleration at the horizontal] 
tail, positive upward 

longitudinal component of velocity of airplane 
positive forward 

wing angle of attack caused by positive (up 
draft) gust velocity and forward speed 

open-loop transfer function, cross feedback 

open-loop transfer function, direct feedback 

closed-loop transfer function 

input toa system 

output of a system 

angular frequency, rad. per sec 

angular frequency, cycles per sec 

differential operator, d/dt 

transfer function, output /input 

longitudinal velocity due to longitudinal force 

longitudinal force due to longitudinal velocity 

longitudinal force due to pitching motion 

longitudinal force due to normal velocity 

longitudinal force due to tail deflection 

longitudinal force due to wing deflection (first 
mode) 

longitudinal force due to wing deflection 
(second mode) 

longitudinal force due to aileron deflection 

longitudinal force due to elevator deflection 

longitudinal force due to elevator tab deflec 
tion 

pitching motion due to pitching moment 

pitching moment due to longitudinal velocity 

pitching moment due to pitching motion 

pitching moment due to normal velocity 

pitching moment due to tail deflection 

pitching moment due to wing deflection (first 
mode) 

pitching moment due to wing deflection (second 
mode) 

pitching moment due to aileron deflection 

pitching moment due to elevator deflection 

pitching moment due to elevator tab deflection 

normal velocity due to normal force 

normal force due to longitudinal velocity 

normal force due to pitching motion 

normal force due to norma! velocity 

normal force due to tail deflection 

normal force due to wing deflection (first mode) 

normal force due to wing deflection (second 
mode) 

normal force due to aileron deflection 

normal force due to elevator deflection 

normal force due to elevator tab deflection 

tail deflection due to tail force 

tail force due to longitudinal velocity 

tail force due to pitching motion 

tail force due to normal velocity 

tail force due to tail deflection 

tail force due to elevator deflection 

tail force due to elevator tab deflection 

wing deflection (first mode) due to generalized 
force on wing (first mode) 

generalized force on wing (first mode) due to 
longitudinal velocity 

generalized force on wing (first mode) due to 


pitching motion 
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= generalized force on wing (first mode) due to Fy, = elevator tab moment due to longitudinal ve 
normal velocity locity 

= generalized force on wing (first mode) due to Fy = elevator tab moment due to pitching motion 
wing deflection (first mode) F's = elevator tab moment due to normal velocity 

= generalized force on wing (first mode) due to M98 = elevator tab moment due to elevator deflection 
wing deflection (second mode) Fy elevator tab moment due to elevator tab de 

= generalized force on wing (first mode) due to flection 
aileron deflection Fir, = elevator moment due to cable force 

= wing deflection (second mode) due to general- For, = elevator tab moment due to cable force 
ized force on wing (second mode) F. = cable force to longitudinal control surfaces 

= generalized force on wing (second mode) due F\;’ = longitudinal force on aircraft due to gust ve 
to longitudinal velocity locity 

= generalized force on wing (second mode) due to ; ag = pitching moment on aircraft due to gust ve 
pitching motion locity 

= generalized force on wing (second mode) due to F;;' = normal force on aircraft due to gust velocity 
normal velocity F,3' = tail force due to gust velocity 

= generalized force on wing (second mode) due F,,’ = generalized force on wing (first mode) due to 
to wing deflection (first mode) gust velocity 

= generalized force on wing (second mode) due Fi’ = generalized force on wing (second mode) due 
to wing deflection (second mode) to gust velocity 

= generalized force on wing (second mode) due to F 33" aileron moment due to gust velocity 
aileron deflection F 53’ = elevator moment due to gust velocity 


aileron deflection due to aileron moment Pn = elevator tab moment due to gust velocity 


aileron moment due to longitudinal velocity q = 6 = angular velocity 


= aileron moment due to pitching motion Ow, = angle of attack of the wing and tail, respec 
= aileron moment due to normal velocity tively 
aileron moment due to wing deflection (first € = angle of wing downwash 
mode) n = tail efficiency = (dynamic pressure at tail 
= aileron moment due to wing deflection (second (dynamic pressure of free stream) 
mode) p = mass density of air 
= aileron moment due to aileron deflection Ww = linear velocity of center of gravity along the 
= elevator deflection due to elev ator moment Te 
= elevator moment due to longitudinal velocity . eo . 
ures ; ° ly = moment of inertia of aircraft about the y-axis 
= elevator moment due to pitching motion : 
: m = mass of airplane 
= elevator moment due to normal velocity f wi 1 tail vel 
. ° . Ses = area of wing and tall, respectively 
= elevator moment due to wing deflection (first wah ; 5 E , : 
lw, lt = distance, measured parallel to x-axis, from ait 


mode) 


= elevator moment due to wing deflection (second plane c.g. to aerodynamic center of wing and 


mode) horizontal tail, respectively, positive rear 


elevator moment due to elevator deflection ward 


= mean aerodynamic chord of the wing 


elevator moment due to elevator tab deflection Cu 
= elevator tab deflection due to elevator tab C,, Cu, Cz = lift, moment, and normal force coefficients, 
moment respectively, of the airplane 
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Buckling of Sandwich Cylinders 
Under Axial Compression’ 


F. K. TEICHMANN,? CHI-TEH WANG,! ann GEORGE GERARD** 
New York University 


SUMMARY 


A theory for the buckling of the sandwich cylinders under axial 
compression is formulated, and the resulting differential equa 
tions are solved. Both the symmetrical and general circumfer 
ential buckling are considered. For long cylinders, the results 
of these two theories agree with each other as in the case of the 
buckling of homogeneous circular cylinders. The buckling 
stress is presented in a nondimensional form in terms of parameters 
that are characteristic of the sandwich cylinders. The theory is 
verified experimentally for sandwich cylinders constructed of a 


core material that is weak in shear. 


INTRODUCTION 


A SANDWICH-TYPE STRUCTURE consists of two ex- 
ternal layers of thin high-strength material and a 
thick internal layer of lightweight material. The former 
are usually called the ‘‘face layers’’ and the latter the 
“core."’ The basic advantage of a sandwich structure 
is its high bending rigidity and consequent buckling 
strength and its relatively light weight. The use of 
sandwich-type structural elements for aircraft construc- 
tion has received considerable attention with the de- 
velopment of high-speed aircraft. At high speeds, the 
surface of the airplane must be smooth without any 
buckling during flight. The use of sandwich construc- 
tion meets this requirement well. 

The buckling and bending cases of sandwich-type 
beams and flat plates have been studied by several in- 
vestigators. The small-deflection bending theory of 
sandwich shells has recently been developed by Reiss- 
ner,' and the case of symmetrical buckling of sandwich 
cylinder under axial compression has been studied by 
Leggett and Hopkins? using the energy method. 

In this paper, the problem of buckling of sandwich 
circular cylinder under axial compression is solved by 

Received May 1, 1950. Revised and received December 19, 
1950 

* The work reported in this paper was carried out under Con 
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obtaining solutions of the differential equations. The 
solution for the buckling load is obtained in a non 
dimensional form that is a function of the parameters 
that characterize the sandwich cylinder. A series of 
tests has been conducted on cylinders that have a core 
material weak in shear. The experimental results are 
found to be in reasonable agreement with those given 


by the theory. 


SYMBOLS 


a = radius of a circular sandwich cylinder 

E = Young’s modulus of elasticity 

G = modulus of rigidity in shear 

h = thickness of the core layer 

k = buckling coefficient of axisymmetric theory 

k. = buckling coefficient of general theory 

l = length of the cylinder 

m = number of half waves along the generator of 
the shell 

M,, My, Mr, = resultant bending and twisting moments in 


the composite shell; unit in moment per 


unit length 


























Fic. 1. Resultant forces and moments acting on a sandwich 
element 
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BUCKLING OF SAN 


n = number of full waves around the circumfer 
ence of the cylinder 


V,, V,, Ney = stress resultants in the middle surface of the 
composite shell; unit in force per unit 
length 

0.0, = shearing stress resultants in the composite shell 
normal to middle surface; unit in force per 
unit length 

t = thickness of a face layer 

x, = coordinates; y is the curvilinear coordinate 


in the direction of the circumference 
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oO. . OU, V. O*w . N, oe “ 
Ox og " (xe a 
(O.M,/Ox) + (OM,,/a00) — Q, = 0 (4) 
(OM,,/0x) + (O0M,/a0d) — Q, = 0 (5) 


where a is the radius of the cylinder. 
The pertinent stress-displacement relationships can 


u,v, U = displacement in x, y, z directions, respectively be obtained from those given by Reissner.' These 
a = (h + t)?/a* relations may be rewritten in the following form: 
3,, B = components of change of slope of normal to 
middle surface of a plate or a shell eins , 4 
7 = (h + t)G.-/Egt | a 2th; k 1 a ( oO " *) 
y = Poisson's ratio ; | — vfLox a O68 a 
d = mra/l iE 
2tE 1 Ov qe 
¢ = N, ‘Fest N ™ f (, Ou L O : ) 
= 2 - é a 
Subscripts | v Ox a oe 
f = core layer . 2tE, 1 Ou Ov 
( = face dayer of the sandwich shell Ny = 2(1 + v,) \a 00 7 or 
» 1 4)2¢ fs / mr.) 
(I} BUCKLING OF CIRCULAR SANDWICH CYLINDERS iw... = 1/26 + HME, (~* _ 1 ) 
UNDER AXIAL COMPRESSION eth. 4 On a OO 6) 
Vv 1/2(h + t)tk, (Og, 03 
Sa : F , aa i 
Equilibrium Equations y = es ae 
By assuming that all resultant forces except N, are 1 2(h + HAE, (108, . 2B, 
extremely small and neglecting the effect of transverse Mi, = 11 4 & , * > 
Z - 4° . . oes os T Vy) ( ¢ ae 
normal stress deformability—1.e., 4, = ©—the equilib- 
rium equations are the same as those used for the homo- Qr = (h + HG 1B, + (Ow Ox)] 
geneous cylinder. The pertinent sandwich element is ry 1 Ow 
i QO, = (h+HG,(8,+—-+ ) 
shown in Fig. 1. ‘ a a 06 
(ON,/Ox) + (ON,,,a00) = 0 (1) et —_ ; : 
Substitution of Eqs. (6) into Eqs. (1) to (5) results 
ON, | ON2, Ov Q, 0 >, in the following equations in terms of u“, 7, w, 3,, and 
aoe Ox * Ox? a ~ 6: 
O°u 1 O°u 1 O*v | Ow 
CO = ie = & = 0 
Ox? a” 06° a Oxod a Ox 
1 O°u 1 O*v _. ov v 1 Ow C7 
(Co + C3) a; = + (C3 + N r) = i = ee oy = . (ap, ) = 0 
a Oxo6 a” Ob" Ox? a” a? O06 a’ 
| Ou 4 in 1 Ov 4 Ly "Ww Cc; O-w Ww 1 O(ag, cz; O(ag, 7 
C2 (¢1 oy —, (C7 Ne. 7 eons ome iy me a + = () i) 
a Ox a” O60 Ox" a* Ob" a’ a@ ON a> O86 
ow C4 0*(aB,) Ce O°(ap, } C7 l 0*(aZ, ) 
— C7 = = : — (aB,) + (C5 + Ce = ( 
Ox a Ox a® O06" a ad, O60x 
v 1 Ow 1 O°(ag,) Cc, O7(aB,) Cs O7(ap, ) C; 
— C7 — > (> cs) — rn T : - ap,) = 0 
a a O6 a’ O6Ox a®— O08 a Ox" a 
where 1 2(h + ¢t)*t Ey — 
Cc = : 0.165 a -+ 4 I: 
¢ = 2E,/(1 — vf) = 2.2ky 7 
Co = 2WtEyy,/(1 — vf) etl 1 20h & DE, 7 
Cs = 2tK,/2(1 + v,) = 0.77 Kz C 1 0.192(h + 1t)*tk 
( T Vy 
1/2(h + t)*tk, i ia 
Cs = 0.55(h + t)7tk, 
l — pv ( (h + t)G 
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Fic. 3. Half wave length of buckling; v = 0.3. 
The numerical coefficients in the above expressions are 
obtained by taking v,; = 0.3. 
Axisymmetric Buckling 
A solution is obtained by assuming 
v= BB, = 0 
u, W, Br = f(x) only 


This is the type of buckling in which the deformations 
are symmetrical with respect to the axis of the cylinder. 
With the above relationships, Eqs. (7) may be reduced 
to one equation in terms of w as follows: 


N.,\ O*w E,t cs\ OPw 
; (: i ) + (¥. i *t) 

C7 Ox! a? C7 Ox? 
Eyt 


a~ 


») 


Z w=0 (3) 
where NV, has been replaced by — N,,. 

For long cylinders, where the boundary conditions at 
the ends are not important, the solution to this homo- 
geneous fourth-order linear differential equation is 


w= > w, sin Ax/a (9) 


m=1 


where \/a = mm 1. 
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Substituting this solution in Eq. (8) gives 


oe) rv! 2 7 2 
y « — Ne = (: +> :) + 
m=1 a‘ a~ 


C7 a 


2E yt C4 Ys 7 AX 
- 1+ - W», sin = () 
a C7 a? a 


For a nontrivial solution, the bracketed coefficient must 
Thus, solving for the critical compressive 


(10) 


vanish. 
loading, 


(11) 
1 + (¢4/c7) (*/a?) 


If a large number of longitudinal waves occur, the value 
of \ which minimizes Eq. (11) is obtained from the con- 
dition 

ON,,/OA\* = 0 
Eq. (12) has two possible real solutions: A? = © 


x2 A+) _ 4 (Et/Ga) + V(1 — v2) 


° ; l (=) (13) 
1 — vf? \Ga 


Note that from Eq. (13) when 
Ept/Ga 2 V1 — vf (14) 


the only real root isA* = ©. 
The normal stress in the faces at buckling is 


Ger = N,,/2l 


c! 


Substituting Eqs. (13) into Eq. (11) and simplifying, we 


obtain 
Cc = k,E,|(h + t)/a} (15) 


where the buckling coefficient , is 


N(h + Z a mn E,t - 
b 4111 — vf)a N(A + 2) 2(1 — v/7)Ga 14 
= | ) 
. E,t (h + t) 
1+ — ie 
2(1 — v/)Ga a 


Values for the buckling coefficient k, and half wave 
length have been computed from Eqs. (16) and (13), 
respectively, by taking vy = 0.38. These data are shown 
in Figs. 2 and 3. When E/t/Ga 2 0.95, the value of 


the buckling coefficient is 


k, = (1/2) (Ga/E,t) (17) 
Under these conditions, the buckling stress 
Oc = (h + t)G,/2t (18) 


has a value that is independent of the wave length. Al- 
though X = © satisfies Eq. (12) under these circum- 
stances, it appears reasonable to consider the wave 
length to be unspecified rather than possessing a van- 


ishing value. 
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General Buckling 

Considering now the general solution of Eqs. (7), 
these equations will be satisfied by assuming the fol- 
lowing relations for the displacements: 


«= } a > Umn Sin nO cos 
1 n=1 


m= 


MAX 


-. MAX 
>) Umn COS NO sin 


m=1 n=1 


x c g . mnrx 
Ww _ i > Wm n sin no sim 1 


m=1 n=1 


(19) 


mmx 
YS Brmn Sin nO cos ; 


m=1 n=1 
MAX 


> Bymn COS nO sin j 


m=1 n=1 


Substituting expressions (19) into Eqs. (17), we obtain 
the following equations: 
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(2.2X? + 0.77n?)tmn + 1.43 mn + 
0.66AW mn = 0 | 
1 ABNNU mn + [2.2n? + (0.77 + g)A2 + ylomn + 
(2.2 + y)NWmn + YABymn = 0 | 
0.66AUmn + (2.2 + y)Ntmn + [Cy + o)rA> + 
yn? + 2.2] mn + YAABrmn + YNABymn = O | 
yAWnn + (0.55adr2 + 0.192an? + y)dBrmn + | 
0.357nAQGBymn = 0 | 


Yma + YNWmn + 0.357nrAAaABrmn + 
(0.55an? + 0.192ad2 + y)dBymn = 0 | 


(20) 


A nontrivial solution will exist only if the determinant 
of the coefficients of these equations vanishes. 

Expanding this determinant, we obtain 
Ag? + Be +C=0 (21) 


where 


A = Nia? (2.2 A? + 0.77n?) }(2) + 0.742 (7) (A> + mn?) + 0.1056 (A? + "| 
a 


a 


B= 


» 


Na? ((*) [1.257(A2 + m2)® + (3.27 A2 + 1.2572) (2 + n2)] + © [(0.192 2 + 0.552) (2.242 + 
a Qa 


0.77n7) + 0.55(A? + m?)? (2.2? + 0.77n?)] + (zy 11.694[(A2 + m?) + 2.62 + n?]} + 0.178 [(A2+ mn?)t + 
a 


(2.62 + n?) (A? + n?)?] + 0.1056 y (2.242 + 0.77n?) [(A2 + n?)® + (A? + n*)'1) 


C=a? (2 [2.52 M4(A2 + m2) ] + 
a 


a 


0.358A4(A? + n*)?] + y[0.1780? 


In general, Ay? may be neglected as compared to the 
Therefore, we have* 


—C/B (22) 


where the negative sign indicates compression. The 
buckling load may be obtained by minimizing Eq. (22) 


other terms. 


a= 


g= 


with respect to \ and n. 

The minimization of Eq. (22) proved to be intractable 
analytically because of the presence of the sandwich 
cylinder parameters a and y in this equation. Conse- 
quently, a numerical treatment of Eq. (21) was at- 
tempted. In the numerical analysis three cases were 
treated: (a) the core considered extremely weak in 
shear with a = 10-4 and y = 2.7 X 107%; (b) the core 
considered moderately weak in shear with a = 104 
2.7 X 10 and (c) the core considered rela- 
* The expansion of the determinant and the numerical calcu 


lations were carried out by R. J. Vaccaro and L. Ting, Graduate 
Students, New York University. 


and y = 


zs 10.932 (A? + n2)4 — 1.863 [(A2 + m2)? — A® + 0.674d2n?] + 


2.2(2.2A? + 0.77n?) (0.192? + 0.5527) — 0.6119A‘n? + 0.4607A2n? — o.os3n+4 + Y) [3.39 A* + 
f 


a 


+ n?)> + (dX? + n?)? (0.465A2 + 0.178? — 0.822A2n? — 


0.357n*) | + 0.1056y7[(2.2A? + 0.77n?)d2(A? + n")'1) 


tively strong in shear with a = 10 
167". 


The values in case (a) were chosen so as to coincide 
with the nominal properties of the cylinders used in the 
experimental study. The other values were chosen as 
convenient multiples of the first set. 


In performing the analysis, all terms were retained in 
the exact expression given by Eq. (21) when they were 
found to be of importance. Thus it was found that the 
contribution of the Ay? term for \ > 10 was to raise the 
value of ¢ by less than a constant value of 4 per cent. 
This was for the most severe condition of m = O and, 
therefore, Ag® can be safely neglected in all cases for 
which \ > 10. 


The results of the numerical analysis are shown in 
Figs. 4, 5, and 6. The values of m and X\ obtained from 


these figures for a minimum ¢ are given in Table 1. 
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In Figs. 4 to 6, curves for only a few values of n are 
shown which are sufficient to ensure a minimum value 
of ¢ for the corresponding \._ It was possible to make 
approximate analytical calculations that served as a 
guide to bracket the minimum m for a particular X. 
This reduced the amount of calculation to the point 


where it was only necessary to calculate from the exact 
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expressions, Eq. (21), the ” curve for a corresponding } 
range. With this procedure, the values of ¢,,;,. are 
reasonably certain, although m and \ may be in slight 
error. 

A simpler approach to the study of Eq. (22) was 
found possible by assuming 


x a f 
ax<xy<!l 


nO) 
( 23) 
It appeared that the conditions given by Eq. (23) would 
cover the region between cases (a) and (b) for which the 
exact solution was obtained. By retaining terms of the 
highest powers in \, the following simplified equation 
was obtained from Eq. (22): 


¢ = — [y(0.932 ” x8 + 0.178 A) + (0.35808 + 
a 


2.52 ™r6§ 4+ 3.39 7 


a a 


1.257 "x® + 1.694 (7) (24) 
Qa a 


Eq. (24) can be put into the following form: 


7) (= 0.178 + 0.9328 ) 
Y y 0.178 + 1.2578 + 1.694? 


; M4) ]/2 [0.178 A + 


oe 


where B = y/ad’. 
To minimize gin Eq. (25), take 


O(¢/y)/OB = 0 
with the result that 


a 0.061 + 0.6038 + 1.579." 
yy? 2(0.178 + 1.2578 + 1.69482)? 


where a/y” = (Eyt/Ga)’. 
Eq. (26) is solved by assuming values for 3 and then 
proceeding to Eq. (25). By substituting the expres- 


sions for g and y, Eq. (25) can be rewritten in the form 


(27) 


Neo = kh + 1G, 


where k, is now the right side of Eq. (25). 

Calculated results for k, and the half wave length oi 
buckling are shown in Figs. 7 and 3. Also shown in 
these figures are the numerical results of the exact solu- 
tion for the three cases discussed previously. From 
the figure, it can be observed that results of the exact 
calculation and the approximation given by Eq. (27) 
are in excellent agreement for all three cases. 

Another interesting feature of the approximate solu- 
tion is that when £ is permitted to be extremely large, 
which presumably corresponds to the case where the 
shear effects of the case are negligible, then Eq. (26) re- 


duces to 
B = 0.813 (y/V @) 
Furthermore, by substituting this value into Eq. (25) 


and simplifying, the critical stress equation with no 
correction for shear effects is obtained 
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o-r = 1.15E,[(h + t)/a] 


The value of the coefficient of 1.15 is only about 10 per 
cent higher than the value of 1.05 obtained from an 
exact solution of the differential equation for this case. 
It should be remarked that such agreement could not 
have been anticipated for the region below case (b) on 
the basis of the conditions given by Eqs. (23). 

The symmetrical buckling theory can be compared 
to the results of the general theory by noting that Eq. 
(15) can be written as 


N.p = 2k, E;t((h + t)/a] 


“ad (28) 
Ner Ft 
= 2k, = k,’ 
(h + t)G, Ga 
where k, is obtained from Fig. 2. The buckling co- 
efficients k, and k,’ are compared in Fig. 7, and it is evi- 
dent that the results of both theories coincide within a 
few per cent. 

An interesting feature of the theories illustrated in 
Fig. 7 is the fact that there are apparently three ranges 
into which the shear effects of the core can be divided: 
(a) When E,t,Ga 2 0.95, the core is extremely weak 
in shear and N,, = (kh + t)G,. (b) When 0.1 < Eyt + 
Ga < 0.95, the core is moderately weak in shear and 
N., is given by Eq. (27). (c) When E,t/Ga < 0.1, the 
core is relatively strong in shear, and the sandwich cyl- 
inder behaves essentially in the same manner as the 
homogeneous cylinder. While shear effects are neg- 
ligible in this region, it appears that nonlinear theory, 
which is necessary to describe buckling of a homogene- 
ous cylinder, may also be required for the sandwich 
evlinder. The linearized solution for the buckling 
stress of the latter can be obtained as 
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o = (1/V 1 — w)E[(h + t)/a] 


The strain energy solution obtained by Leggett and 
Hopkins?’ also agrees with the results shown in Fig. 7 
when their solution is simplified for a core that does 
not carry any axial load. A slight difficulty in their 
solution which was not noted by Leggett and Hopkins 
is that, in the region A—A’ shown in Fig. 2, the buckling 
stress given by their equation goes to zero, for in this 
region \ becomes imaginary for this branch and, con- 
sequently, the \ = © branch must be used. 


(II) EXPERIMENTAL RESULTS* 


Coincident with the theoretical investigations given 
in the previous sections, an experimental program on 
sandwich cylinders under compressive end load was 
conducted to provide data for correlation with the theo- 
retical results. Acquisition of experimental data by 
itself was considered a worth-while objective, since test 
information on sandwich cylinders is completely lacking 
in the literature. In addition, these data prove help- 
ful in providing an inside into the 
of the effect of various terms encountered in the theo- 


‘ 


‘order of magnitude” 


retical development and thus permitted considerable 


simplification in certain instances. 


Test Program 


In an effort to test cylinders within a range of prac- 
tical interest, it is desirable that the radius to thickness 
ratio should be as large as possible. Whereas a rela- 
tively large value of the a/t parameter is easily obtain- 
able for a homogeneous cylinder by constructing the 
cylinder of sheet of several thousandths of an inch 
thickness, the sandwich cylinder presents the problem 
that considerably greater thickness is required for the 
sandwich for practical reasons. 

A '/s-in. thick sandwich was considered the mini- 
mum thickness feasible, and consequently rather large 
radii were chosen so that the a/f values would be within 
a range of possible use in aircraft construction. Three 
different radii were used to permit a variation in the 
a/t parameter to be studied within the capacity of the 
testing facilities. To determine the length effect, if 
any, two different lengths were used in the investiga- 
tions. 

All other details of the specimens were similar. The 
24S-T aluminum-alloy faces of 0.010-in. thickness were 
bonded to a !/s-in. cellular cellulose acetate (CCA) core 
of approximately 4.5 to 5.0 Ibs. per cu.ft. density. At 

* Acknowledgment is made to R. Freeman, H. Hilton, and J 
Tolciss, Graduate Assistants in the Department of Aeronautical 


Engineering, for participation in the experimental program 


TABLE | 
a ¥ Ya) n r Cmin./¥ Va/y 
10 2.74 x 10 2.7 x 0 x 1.00 3 7 
10-* 3.4 xX wr? |. xX 10 0 18 0.604 0 37 
10~4 2.7 x we 2.04 X 10 7 7 0.0755 0.037 
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TABLE 2 
Presentation of Test Data 
Calculated 
Nominal Inside Total Face Core Ultimate Drop 
Length, Radius, Thickness, Thickness, Thickness, Load, Load, 
Cylinder * In. In. In. In In. Lbs. Lbs 
1012At 12 10.05 0.144 0.0107 0.1226 9,675 4,225 
1012B 12 9.99 0.134 0.0107 0.1126 11,825 4,025 
1024A 24 10.09 0.138 0.0107 0.1166 10,625 
1024B 24 10.05 0.135 0.0107 0.1136 13,985 
1212A 12 12.00 0.142 0.0107 0. 1206 18,175 
1212B 12 11.89 0.136 0.0107 0.1146 18,475 8,950 
1224A 24 11.88 0.1382 0.0107 0.1106 17,625 8,225 
1224B 24 Premature failure due to poor bonding 
1412A 12 13.89 0.143 0.0107 0.1216 26,950 2,250 
1412B 12 13.86 v. 134 0.0107 0.1126 33,800 2,400 
1424A 24 13.88 0.132 0.0107 0.1106 21,150 5,900 
1424B 24 13.84 0.131 0.0107 0.1096 18,450 7,100 


* Note that the cylinders are identified by the radius that is given by the first two digits and the length that is given by the last 


two digits. 


+ This cylinder was initially buckled, during manufacture, at the location of buckling under load 


each end of the cylinder, a 1'/:-in. deep hardwood in- 
sert replaced the core to allow for the uniform introduc- 
tion of end loads into the cylinder. 

The cylinders were manufactured by Skydyne, Inc., 
Port Jervis, N.Y., and their fabrication details are 
quoted: ‘‘All aluminum faces are coated with a metal 
primer prior to assembly to the core. The primer is 
an adhesive manufactured by the Casein Company of 
America and bears their code number of NT-442. 
The CCA core is then planed to the thickness of !/s in., 
and solid wood reinforcements are glued to the edge. 
A secondary adhesive is used to combine the aluminum 
faces and the core; the adhesive used is Plastron 250-2. 
Although this resin adhesive will polymerize at a room 
temperature of 75°F., sufficient heat was applied dur- 
ing the fabrication of the cylinders to accelerate the 
‘setting’ of the bond. Pressure necessary to bring the 
surfaces in intimate contact was applied, and the as- 
sembly was molded in a period of 2 to 3 hours. The 
final operation consists of trimming or finishing the 
ends of the cylinder to assure an even plane and square- 
ness to the ends.”’ 

Load was applied to each specimen by means of a 
Baldwin-Southwark universal hydraulic testing ma- 
chine. The 10-in. radius cylinders were tested in a 
machine of 60,000-lb. capacity, whereas the remainder 
were tested in a 200,000-lb. capacity machine. The 
specimen rested on a '/4-in. machined plate, which in 
turn rested directly on the base of the machine. A 
1'/»-in. machined steel plate was centered on top of the 
cylinder and load was applied to this through a self- 
aligning compression head. 

During the progress of the test, load was applied 
slowly. At appropriate increments of load, strain-gage 
readings at four stations and deflection readings at 
three stations were taken. Near the buckling load of 
the cylinder, the load increments were reduced. At 
failure, both the maximum load and the drop load at 
which the machine stabilized were recorded. 


Ultimate Loads and Dimensions 


A summary of the important test results is given in 
Table 2. 
are those that were obtained prior to test at the esti- 
mated center of the buckle pattern. It was noted, how- 
ever, that this location usually did not coincide with the 
location of minimum thickness or maximum length of the 


The various dimensions given in this table 


cylinder. 

Material property tests conducted on the CCA core 
material indicated that a shear modulus of G, = 2,000 
Ibs. per sq.in. + 20 per cent was a representative value. 
Since the stress in the faces at buckling was well below 
the proportional limit of the material, no material 
property tests were conducted on the face material. 


Description of Failure 


In general, the mode of instability observed on the 
cylinders tested can be roughly divided into three 
groups: (a) Cylinders in which buckles appeared in 
isolated locations before snap buckling precipitated 
failure: 1012A, 1012B, 1212B. (b) Cylinders in which 
no visible evidence of buckling occurred before snap 
buckling precipitated failure: 1024A, 1024B, 12124, 
1224B. (c) Cylinders in which outward bulging either 
preceded or accompanied snap buckling: 1412A, 1412B, 
1424A, 1424B. 

The cylinders in group (a) exhibited isolated buckles 
that grew as the load was increased. These buckles 
usually formed the nucleus of the snap buckle pattern 
that occurred at failure with explosive violence. At 
failure, the ultimate load dropped to a value approxi- 
mately 40 per cent of the peak value. The behavior of 
group (b) was similar with respect to the drop in load. 
The characteristic buckle pattern obtained in these 
groups is shown in Fig. 8. 

Group (c) cylinders exhibited a short wave length 
outward bulge, usually at the junction of the CCA core 
and hardwood insert, in addition to a snap buckling 


pattern. The former mode of instability appeared 
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Typical cylinder after test 


to be characteristic of the type associated with wrinkling 
of the faces. In the cases where the outward bulging 
preceded failure, it was noticed that portions of the bulge 
area changed to the characteristic snap buckle pattern 
at, or immediately after, failure. In addition, the 
drop load was much less than for the other two groups, 
averaging about 20 per cent of the peak value. 


Buckle Patterns 


Data on the buckle patterns obtained in the tests are 
given for each cylinder in Table 3. 

The shape of the buckles obtained from the snap 
buckle patterns were diamond-shaped as shown in Fig. 
8 and of approximately 3 to 1 aspect ratio as given in 
Table 3. They indicated the characteristic tendency to 
buckle inwardly. It is noted for purposes of compari- 
son that, for a homogeneous cylinder, the buckle is 
usually diamond-shaped and of an aspect ratio close 
to unity. 

Another detail of importance was the occurrence of 
only one buckle in the axial direction for all the cylin- 
ders. In testing Cylinders 1012B and 1024A, excessive 
deformation was applied to the cylinder to determine if 


TABLE 3 
Data on Buckle Patterns 


Buckle Size 
Circum 
ferentially 


Longitu- Number of 


dinally eae Circum- 

L, = I/m, i n ferential Lo _» 
Cylinder in in. Buckles *  & n 
1012A 1.5 4.0 24 2.67 
1012B 1.5 4.0 24 2.67 
1024A 1.5 1.0 24 2.67 
1212A 1.5 1.5 24 3.00 
1224A 1.5 4.5 24 3.00 


* Note that these buckles were sinusoidally arranged around 
the circumference (see Fig. 8) 
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additional buckles would form in the axial direction. 
It was observed that no further buckling occurred. 
This behavior is probably due to the formation of a 
weak area at the buckle location once buckling occurs. 
Further deformation of the cylinder appears to be com 
pletely restricted to the buckle location with ultimate 
crushing of the core and separation between the face and 


core. 


Influence of Length 

A 12- and 24-in. length of each cylinder were included 
in the test program to determine the influence of length 
upon the buckling characteristics within this range of 
length. Allowing for the hardwood inserts in the end 
of each specimen, the free length, L,, of the cylinders in 
which buckling is likely to occur is reduced to 9 and 21 
in., respectively. 

The buckling load per inch, N,, as a function of 
L,/a, is shown in Fig. 9. Within this range of the 
parameter and considering the scatter of the test data, 
no conclusions as to the effect of length could be deter 
mined. Examination of the buckle dimension data in- 
dicates, however, that no effect of length on the buckling 
load should occur, since, for the shortest cylinder, the 
axial wave length was approximately one-sixth of the 
free length. It is anticipated that only when the free 
length is several times the axial wave length will the 
length influence the buckling characteristics. 


Physical and Material Property Variation 


Examination of the dimensional properties of each 
cylinder indicated that the variation in dimensions from 
the average were of the following approximate orders of 
magnitude: (a) face thickness, ¢ = 0.0015 in. (com- 
mercial tolerance); (b) core thickness, + 0.008 in.; 
(c) inside radius; a + 0.030 in.; (d) length, 2 + 0.008 
in. 

Variation in the shear modulus of the CCA core was 
found to be G, = 20 per cent. The variation of /, of 
the face material is known to be negligible for aluminum 
alloys. <A value of Fk, = 10.5 X 10° Ibs. per sq.in. was 
used. 

An attempt was made to correlate the location esti- 
mated to be responsible for buckling with initial imper- 
fections such as minimum values of thickness or maxi 
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= | dD It appears that, for this range of E,t/G.a, the critica] 
| | stress is independent of the radius a and depends essen- 

300 meee pai | tially upon the core properties. 
x, W,etnene, rt d To test the validity of this conclusion, the test data 
IW8/i0 oof 2 jo ‘|6 «1600 |psi nA are plotted in Fig. 10 asa function of a. Also shown are 
0, the results obtained from Eq. (29) based on the upper 
[ and lower limits of the shear modulus, G,, given pre- 
= SS Se viously. Since no correlation of failure location with 
dimensional imperfections was obtained as discussed 
‘ in the previous section, it is assumed that the variation 
° ® 0 " 2 % ” '§ in the shear modulus is primarily responsible for buck- 
a” ling. 

Fic. 10. Comparison between test data and theory based on While the test data in Fig. 10 appear to indicate an 


limits of variation of G,. 


mum values of length. From the tabulation of data 
given for each cylinder in reference 6, no correlation was 
indicated. Since correlation with dimensional imper- 
fections was lacking, it is assumed that the buckle loca- 
tion is caused primarily by local weak spots created by 
the variations in the core shear modulus. The fact that 
only one longitudinal buckle formed in all cylinders 
appears to be further confirmation of this hypothesis. 
Additional evidence of this type of behavior of sand- 
wich elements subject to buckling has previously been 
offered by Hoff and Mautner.* They suggest that, in a 
sandwich structure subject to uniform load, such as a 
beam in bending, the results are a statistical average 
of the core properties. On the other hand, a structure 
subject to buckling responds to local weak spots caused 
by defects in the core properties, and, hence, the results 
should be based on the lowest properties rather than on 


the average properties. 


Correlation of Theory and Test Data 


It was found herein that for values of Kyt/G.a 2 0.95, 


the critical loading is given by the following equation: 


N,, = (hk + 0G, (29) 


apparent influence of the radius, it is felt that this ap 
parent effect can be explained by the variation in the 
shear modulus. Consideration must also be given to 
the fact that the 10-in. radius cylinders were tested in a 
60,000-lb. capacity testing machine, whereas the re- 
mainder were tested in a 200,000-lb. machine. This 
may have some bearing on the lower values of the 10-in. 
radius group. Furthermore, the cylinders were manu 
factured in order of increasing radius, and, if some allow- 
ance is made for an improvement in manufacturing 
technique with experience, then the large radius cylin- 
ders are favored. 

Considering all these factors, it appears that reason- 
able correlation between the critical load obtained 
from the theory and the test data was obtained within 
the region of /:,t/G.a of the experiments. 


REFERENCES 


' Reissner, E., Small Bending and Stretching of Sandwich-T ype 
Shells, N.A.C.A. T.N. No. 1832, March, 1949 

2 Leggett, D. M. A., and Hopkins, H. G., Sandwich Panels and 
Cylinders Under Compressive End Loads, British Aeronautical 
Research Council Reports and Memoranda, No. 2262, 1949 

3’ Hoff, N. J., and Mautner, S. E., Bending and Buckling of 
Sandwich Beams, Jounal of the Aeronautical Sciences, Vol. 15, 
No. 12, pp. 707-720, December, 1948. 











dui 
Th 
str 
str 
an 
col 
ma 
mi 


of | 


Du 
sid 


al 


pre 


to 
tai 


me 


tia 


to 


ritical 


-ssen- 


data 
Nl are 
Ipper 
pre- 
with 
ussed 
ation 
uCck- 


fe an 
S ap 
1 the 
n to 
ina 

re- 
This 
Q-in. 
anu 
low- 
ring 
ylin- 


Son- 
ined 
thin 


ly pe 


and 
tical 


g of 











A Minimum Energy Solution and an 
Electrical Analogy for the Stress 
Distribution in Stiffened Shells’ 


L. A. GORANt 
McDonnell Aircraft Corporation 


SUMMARY 


The governing equations that can be used in numerical proce 
dures for the stress analysis of stiffened shells are developed. 
The development is based on the usual assumptions of airplane 
stress analysis—that the ribs are rigid in their own planes, that 
stringers carry only axial load, that webs carry only shear load, 
and that average shear flows between ribs are obtainable from a 
consideration of axial-ldad differences in the stringers. The 
mathematics of the solution, which is based on the principle of 
minimum strain energy, is facilitated by the use of the method 
of Lagrangian multipliers 


The problem considered is a generalization of that solved by 
Duberg,! in that effects of taper and curved shear webs are con 
sidered, while Duberg’s treatment excludes the effects of taper 


and considers only flat webs. 


A simple electrical analogy to the stiffened shell problem is 
presented, and the method for its application is briefly discussed. 


A numerical example illustrating the application of the method 
toa swept plate is presented in the Appendix. The solution is ob- 
tained by the use of the governing equations and also by a nu- 


merical solution of the equivalent electrical circuit. 


SYMBOLS 


A = area enclosed by web median line and lines 
joining the ends of the web to the torsional 
reference axis (see Fig. 3), sq.in. 


A = cross-sectional area of stringer, sq.in. 

E = modulus of elasticity, lbs. per sq.in. 

F, = net external load along y axis at rib station 
7; positive load produces tension, Ibs 

G = shear modulus, lbs. per sq.in 

I = longitudinal current, amp. 

1 = transverse current, amp. 

1 = rib station when used as a subscript 

j = subscript denoting stringer or web element 

L = length of stringer or web element, in. 

M = external bending moment about x axis at 


rib station 7; positive moment produces 
compression in the upper surface, in.Ibs. 

M = external bending moment about zs axis at 
rib station i; positive moment produces 


compression in trailing edge, in.lbs 


P = stringer axial load, lbs. 

q = web shear flow, Ibs. per in 

R = longitudinal resistance, ohms 
r = transverse resistance, ohms 
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li = external torque about torsional reference 
axis at station 7 — 1; positive when clock 
wise, looking toward the root of the struc 


ture, in.Ibs. 


t = sheet thickness, in. 
U = strain energy, in.lbs. 
Si = function formed in applying Lagrangian 


multiplier method, in.Ibs. 


e = voltage, volts 
w = developed width of shear web element, in. 
=, % % = cartesian coordinates (see Figs. 1 and 3), in 


X, Y, Z,@, A = Lagrangian multipliers 


INTRODUCTION 


thee CUSTOMARY METHOD for determining internal 
loads in the stiffened shells employed in aircraft 
structures is to use ‘bending theory,’’ which assumes 
that transverse sections originally plane before bending 
remain plane after bending. It has long been recog- 
nized that this method of analysis gives close to correct 
results only when the effects of shear distortion are 
small and that there are a number of frequently occur- 
ring cases where the use of the method gives appreciable 
errors. 


Problems in which the shear deformation effects are 
considered are usually called ‘‘shear lag’’ problems, and 
a number of solutions are available for such typical 
separate problems as the stress distribution near the 
root of a four-element box beam subjected to torsion, 
the stress distribution in shells in the region of a cut- 
out, and the stress distribution in flat panels subjected 
to local concentrated loads. 


Since, in practice, many of these problems occur 
simultaneously in the same structure, it is desirable to 
have a method that can handle all of these ‘‘shear lag”’ 
effects in the same solution. One important type of 
structure where such a method is a necessity if correct 
loads are to be obtained is the sweptback wing. 


Duberg,' using basically a method of consistent de- 
flections, has developed a numerical method for un- 
tapered structures having flat webs which handles all 
of the shear lag effects in one solution. In the present 
paper, the governing equations for solving the same 
problem, including the effects of taper and curved 
presented. The compatibility 


webs, are necessary 
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equations are obtained by means of an energy method 
combined with the method of Lagrangian multipliers. 

The method here employed for deriving the com- 
patibility relationships has a number of advantages over 
the methods usually used in solving problems in inde- 
terminate structures. The use of an energy method 
eliminates the often difficult task of correctly account- 
ing for all deformations which is required in a consistent 
deflection method. The use of Lagrangian multipliers 
simplifies considerably the algebraic manipulations 
that are required in the usual minimum energy analysis. 
As illustrated in the example presented in the Ap- 
pendix, the method can handle peculiar boundary 
conditions, such as are found at the roots of sweptback 
wings. 

In order to obtain the internal load distribution in 
the usual stiffened shell structure by the method pre- 
sented in this paper, it is necessary to solve a large num- 
ber of simultaneous equations. An iteration procedure 
is probably the only practical means of accomplishing 
this by analytical methods. The use of the electrical 
analogy presented herein would provide a rapid and 
direct solution to the problem. 


DERIVATION OF EQUATIONS 


Basic Assumptions 


In accordance with usual airplane stress-analysis 
practice, the actual structure is replaced by an idealized 
one equivalent to that shown in Fig. 1. The analysis 
is based on the following assumptions: (1) Stringers 


carry only axial load; (2) webs carry only shearing 
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Fic. 2. Equilibrium of stringer element 


force; 
panel”’ effects are negligible; (4) ribs are rigid in their 
own plane but have zero transverse stiffness; (5) taper 
‘moderate’ and is such that the 


(3) the geometry is such that ‘‘skewed shear 


of the structure is 
stringer centerlines all meet at a common point; (6) 
Hooke’s law applies. 


Numbering System and Sign Conventions 


The numbering system and sign conventions em- 
ployed are shown in Fig. 1. Rib stations are assigned 
even numbers, with the root station being number 0. 
The stringers are numbered consecutively starting 
with number one. A tensile loading is considered 
positive. 

Looking inboard, positive shear flows are assumed to 
act clockwise for a left-hand beam, such as shown in 
Fig. 1. The numbering system for shear flows is de- 
scribed in Figs. | and 2. 


Method 


The procedure used for the determination of loads 
throughout the structure is to write the set of equations 
requiring that the strain energy stored in the system 
be a minimum. In addition to satisfying the require- 
ment of minimum energy, at each station the require- 
ments of statics must be satisfied. The mathematical 
problem of finding the loads may be described as a 
problem of a constrained minimum, with the con 
straints being the conditions of statics. In solving 
this problem the Lagrangian multiplier method? is em- 
ployed. 

The basic formulas will be developed in the region of 
stringer j at rib station 7. 
Equations of Statics 

Stringer Equtlibrium.—Since the internal stringer 
forces must equilibrate the external loads, at each 


rib station the following three relationships must 
hold: 
d? = > (P,,,) — F; = 0 (1) 
J 
gd) = > (x, ,P;,.) + Mz: = 0 (2) 
} 
o3°) = >> (8): P),1) + Me, = 0 (3) 


) 


where the summations include all of the stringers. The 
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STRESS DISTRIBUTION 
symbol @ is here introduced for convenience in the sub- 
sequent analysis. 

Shear Flow Relationships.—In writing equations for 
shear flows, the “‘unit’’ method of beam analysis* is 
used. From Fig. 2, Eq. (4c) is obtained. By suitably 
changing subscripts, Eqs. (4a), (4b), and (4d) are then 


obtained. 
o4'? i. % = Li (qj, 3 ; = qj bw 1) + 
Pya,i — Py,i2 = 0 (4a) 
os’ ie = Lislqy, in — Qj 1141) + 
Pya,ize — Pja,i = 0 (4b) 
gs” ' = Lilgj4i, i-1 — Gj, i) + 
P;,i — Pj,i2=0 (4c) 
yb? *T? = Loan (G541, i41 = qi, i+1) + 


Pi we —P;,=0 (4d) 


Torsional Equilibrium.—tIn general, for a tapered 
beam, in writing the equations of torsional equilibrium, 
it is necessary to include the effect of the components 
of stringer load which are coplanar with the shear 
flows. In the present analysis, as shown in Fig. 3, it 
is assumed that the stringers all meet at a common 
point. The reference axis for torsional equilibrium is 
taken through the projection of this point, and, con- 
sequently, the stringer load components need not be 
considered. 
For bay 7 — I, 


d;‘' ; = >> (q;, +-1°2A);, 1) — Fa = (Sa) 


and, for bay7 + 1, 


ost) = > (gy, i41°2Ay, in) — Tins = 0 (Sb) 
U 


where A,, ;; is the area, in the plane of the paper, en- 
closed by the median line of web j, 7 — 1 and lines 
joining the ends of the web to the torsional reference 
axis and where 7;_; is the external torque about the 
torsional reference line. Clockwise torque looking 
inboard is considered positive. The summation in- 
cludes all of the web elements in the bay. 


Energy in Structure 
The total elastic energy stored in the structure is com- 
posed of the cap (stringer) axial energy plus the web 


shear energy. 
U = Ucaps + Uwers (6) 
For cap j, between stations 7 — | and7z + 1, the ex- 
pression for energy is taken to be 
P,, 7L;/2A;, iE 


; is the cap area, and 


where L; = (Lin + Lit1)/2, A; 


Eis the modulus of elasticity. 
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Station 7 — 1. 


Fic. 3. View looking inboard 


For the entire structure the expression for cap energy 


becomes 
Cues ™ pM 2 (P;, 7L:/2A;, EF (7) 


and the summation includes every cap element. 


Similarly, the expression for web energy for web j, 


Gi. 4-2" (“) Lin 
2 t/;,+i G 


where w is the developed width of the web, ¢ is the web 


i — lis 


thickness, and G is the shear modulus. 


For the entire structure the web energy is 
- gj, i-1" Ww a 1 
Um = TEeer (2) sw 
. 2 on t j, i-1 G 


and the summation includes all web elements. 
Substituting the results of Eqs. (7) and (8S) into Eq. 
(6) gives 
Pe tig ( 2 /w L 
U=TE s+ ret=(%) = © 
co. 2A j, iE i? 2 t 1 G 
Minimum Energy Equation 


In comparing the preceding equations with the dis- 
cussion in reference 2, it is seen that Eq. (9) gives the 
function that is to be minimized, and Eqs. (1) through 
(5) express the necessary constraints. 


The function L’* is then formed as 
U*=U+ DY + DX igo? + 
> Zibs? + Do Ay, hi? + VOu@s''-" (10) 


where the quantities, }’;,, X;, Z;, A;,;, and 9, are the 
Lagrangian multipliers that correspond to Eqs. (1) 
through (5). 

The equations for minimum energy are obtained by 


writing 


oU* OP;,; = 0, 0U*/dq;. i1 = 0 (11) 
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Writing Eq. (11) for stringer loads P;,,; and P,;, ; and for web loads q,, ;-1 and qj, :4: gives — 
rey GF ’ Od, ‘” Od»‘” O¢;°” Og4'/ 1s 0¢,°/ 1, i+2) 
+ Y;, ».¢ Li oe A; aS + A —l, i+2 = 0 
oP; 1, i OP; -l,i * oP, ., 7 sa OP, l,i OP; ‘¢ - OP; 1, i | 
(12 
aU Og Ogu! Og" Ogi! 9 Odal* | | 
} i xX i Sis Zi A;, i = A i+2 = 0 
a. ae a Oe, oP, 
OU ra) (j-—1, i ra) (j, i fa) 5! 1) \ 
ttuce +a thee = 0 
Og), i—l Og), i—l Og), i—l Oq), i—1 
(13) 
OU Dd, (F-b §+2) Od,(s +2 Od; (i+! | 
+ A; 1, i+2 ae + Aj, i+2 ~ + Oi41 - = | 
Og), i+1 Oq;, i+l° Og j, i+l1 Og), i+1 


Typical terms of Eqs. (12) and (13) are evaluated below: 


QU/OP 1, ; = Pia,; (L;/Aj-1.:E), 


From Eq. (7), 


OU/OP,.; = Pi,« (Lo As iB) 








ou (“) rey Gi (“) List (14a) 
= qj, i-1 “et = qj, i+1 . 
Og), i 1 fs t/ 5, 1-1 G Og), i41 Ls 3. 441 G 
From Eqs. (1) through (5) 
0g," nee Od."' acd O¢;"' me 
oP, ; a a 
ra) cz. $) Og,‘ * +”? Og,’ l, 2 
ae = -1 «= hes (14b) 
oP. i oP; i Od), i—l 
Od (J, t) Od; ' i 1 a 
= Li ‘5 _ 2: 1 j, i-l 
Og), i—l Og), = 
Putting the results of Eqs. (14a) and (14b) into Eqs. (12) and (13) gives is 1 
I eqt 
Pii gt Vit tai Xi t+ ahi + Apa — Maine = 0 (15a) 
44 3-1, 12 
L; . , 
| { E + } i + Xj, i X + 3), iZ i + A;, tae Aj. iw. = 0 (15b or 
Ww Bes 1 “- 
qj, i-1 — bt Aga, in — Ay: Lin + O-1°2A;, 4 = 0 (16a 
iiaa & 
w Lis 7 hei 
qi, #1 SH Ay, te Linn — Ay, tg2 Linn + Oi41°2A;, in: = 0 (16b) . 
t/ si G tio 
By combining Eqs. (15a), (15b), (16a), and (16b) the Lagrangian multipliers Y;, Aj—1, ;, 4j-1, i42, Aj, i, and = 
Aj, i42 may be eliminated; this leads to the final equation to 
iL w l >. Oi i Sig 
Pr i +X i X; + 3; iZ i = a qj, i-1 ( ) ( + 2A i—1 sel P; ‘4 “ie “18 
ae. wet. ant ae 
x Z (“) (“.) 7 oe ae 
Xx pet — Sy, ee GG, 448 ‘. SFA , i+ = (14 
= d tf3. pa NG Lis 1 
an tr 
At a cutout or a built-in section, Eq. (17) does not probably the only practical means of obtaining a solu- Str 
apply. However, the modified compatibility relation- tion by analytical methods. ‘ii 
ship required in these regions may be obtained by a An example illustrating the application of the method 
- . if ‘s - ‘- ° ° 2 She 
procedure similar to that used to derive Eq. (17). to the determination of loads in a swept plate by means 
Sets of Eq. (17), or its equivalent, together with sets of an iteration procedure is given in the Appendix. aa 
- . - . ~ - . 
of the equations of statics, Eqs. (1) through (5), are 
then sufficient to- determine the loads throughout the ‘Si#nificance of the Lagrangian Multipliers Str 
structure. Because of the large number of equations It is of interest to consider the physical meaning of 
needed for the usual structure, an iteration procedure is the Lagrangian multipliers 6;: and 6;,:. If Eq. (16) ate 
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Fic. 4. 


is written for every web in bay 7 — | and this set of 
equations is added, the resulting equation is 


+ 01 >2A;,i1 = 0 


> | (“) | 
j wa - t j, i-1 G } 


or 


hence, the quantity —@;; represents the relative rota- 
tion of rib 7 with respect to rib 7 — 2. Similarly, 
—6,,; represents the rotation of rib 7 + 2 with respect 
to rib 7. 

For the multipliers X; and Z;, no such simple physical 


significance has been found. 


TABLE | 
Correspondence Between Structure and Electrical Network 


Structure Network 
Item Symbol Item Symbol 
Stringer load Pj, i Longitudinal current Tj, i 
; L; 
Stringer flexibility Longitudinal resistance Rj, i 
Aj, ik 
Shear load qj, iaLi-a Transverse current tj 
u l : 
Shear flexibility - Transverse resistance rj 
tJj,i-1 LiaG 
Lagrangian xX; Voltage factors XxX 
multipliers £ Zi 
6-1 6 
Stringer coordinates xj, ; Longitudinal voltages Xj, i 
3j. i Zi, 6 
: 2A;, ’ 2Aj, i 
Web geometry term lransverse voltage P 
I “i 
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ELECTRICAL ANALOGY OF AA PORTION OF STRUCTURE 


Correspondence between physical structure and equivalent ‘‘direct current’’ electrical network 


ELECTRICAL ANALOGY 


Direct-Current Analogy 


Eq. (17) indicates a simple correspondence between 
the physical structure and an equivalent electrical net- 
work. This network is shown in Fig. 4. Since the 
circuit is purely resistive, its voltages could be either 
direct or alternating; however, for convenience, this 
circuit is called the ‘“‘direct-current”’ circuit. 

The term-by-term 
physical structure and the electrical network is given in 
Table 1. 

As shown in Fig. 4, in each longitudinal arm, two 
external voltages x;, ;X; and z;,,Z; are introduced. In 


correspondence between the 


each transverse arm, the voltage (24;, ;-1/L;-1)6;-1 is 
Across each bay, the voltage factors 


These factors, which de- 


introduced. 
X,, Z;, and 6;—; are constant. 
pend only on the external loading on the structure, 
must be adjusted (by trial and error) until the currents 
are such that the equations of statics [Eqs. (2), (3), and 
(5)] are satisfied. The resulting currents throughout 
the network then give the loads throughout the struc- 
ture. 

It is of interest to note that the electrical analogy 
presented by Newton‘ is a special case of this analogy. 
The structural problem considered by Newton is that 
of a flat symmetrically loaded panel. For this case 
the Lagrangian multipliers are zero throughout the 
structure, and the network presented herein becomes 
identical to that presented by Newton. 
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Alternating-Current Analogy 


The alternating-current circuit shown in Fig. 5 pro- 
vides a direct method for obtaining the loads, starting 
with known external moments. This circuit is identi- 
cal to the direct-current circuit except that the volt- 
ages are obtained from the secondaries of a number of 
transformers. The turns ratio of each of the trans- 
formers is proportional to the physical geometry of the 
The transformers used must be essentially 
perfect. Losses and phase shifts must be negligible in 
order to obtain an accurate analogy. 

Fig. 6 shows typical transformer circuits for the 
transformers associated with terms 6;_; and X,. Con- 
sider, for example, the X; circuit. The primary volt- 
age (initially unknown) is X;. 
stringer j, 1, the secondary voltage is x;, ;X;. 
current /;, ;, in the secondary of this transformer, the 
Since the primaries of the 


structure. 


Corresponding to 
For a 


primary current is x;, ;J;, ;. 
transformers are in parallel, the total input current 


1S 


Es: 1 = > x;, lj, i 


which, from Eq. (2), is seen to be —.V/, ;, the bending 
moment acting on the section. Hence, by reading the 
current /, ;, the proper voltage X; can be impressed 
such that the internal loads equilibrate the external 


loads. 


Discussion 


The use of the direct-current analogy would reduce 
the number of unknowns to be determined from 2n 
per bay (where 7 is the number of stringers) to three per 
bay (X;, Z;, 0::). The alternating-current analogy 
provides a direct means for obtaining loads; however, 
its use may prove to be impractical because of the 
high cost of ‘“‘perfect’’ transformers. 
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Fic. 5. Secondary circuit of alternating current analogy. 
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Typical primary transformer circuits 


CONCLUSIONS 


The use of the principle of minimum strain energy, 
together with the method of Lagrangian multipliers, 
provides a convenient method for determining the loads 
in indeterminate structures in which consistent deflec- 
tion methods are difficult to apply and it is desired to 
avoid the unwieldy mathematics associated with the 
usual minimum energy method. In this paper, the 
method is applied to the problem of stress distribution 
in stiffened shell structures, and the general compati- 
bility equation needed for the solution of this problem 
is obtained. As illustrated by the example presented 
in the Appendix, the method handles complex boundary 
conditions such as those that arise in the stress analysis 
of sweptback wings. 

The use of an electrical-network computer, based on 
the analogy presented, should provide a practical and 
rapid means for determining loads in complex stiffened 
shell structures. 


Appendix 


NUMERICAL EXAMPLE 


As an example, the load distribution in a skewed 
stiffened plate subjected to a uniform shear flow of 100 
Ibs. per inch on its boundaries, as shown in Fig. Al, is 
determined. This problem occurs in the stress analysis 
of a sweptback box beam having only end ribs. The 
loads are obtained by two methods: (1) direct solution 
based on Eq. (17) and its equivalents, and (2) solution 
by numerically determining the currents in the equiva- 
lent direct-current electrical circuit. 


Direct Solution 


Assumptions.—In analyzing the structure, a series 


of equally spaced transverse stiffeners (hereafter called 
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ribs), as shown in Fig. Al, are assumed to exist. As 
noted in the derivation of the general compatibility 
equation, these ribs are assumed to have infinite axial 
stiffness. In addition, the end skewed ribs are as- 
sumed to have infinite axial stiffness, and all ribs are 
assumed to have zero transverse stiffness. The end 
rib is assumed to transfer its load to the pane] as shown 
in Fig. A2. 


Ro,2 + Rs.4 + Rao + Rss = 40/2 


Hence, from statics 


where the loads R are measured in kips. 
In terms of the cap loads, this gives 


Po » + P34 + Ps. 6 + Ps, 8 = 4.00 (Al) 


Analysis Procedure.—For a flat panel, the Lagrangian 
multipliers Z and 6 do not appear in the compatibility 


relationship, and Eq. (17) becomes 


P cL 4 (*.) P ‘. 
it 4 E qi, i-1 1G Rs j-l,i A; Po 


jy 
u 


qi, i+) ( ) + X j(x;, j — xX} 1, i) = () (A2) 
{G j, i+1 


For an iteration procedure, the compatibility relation- 
ship can be expressed in a more convenient form. 
Writing the equivalent of Eq. (A2) for the next web 


gives 
| 4 (“.) 
Qi+, i , “ 
a ee 


Pj, i ” = > Qj+i, i+1 (*.) 
: A;. iE tG j+i, i4+1 4 


X (xj41,14 — Xj.) = 0 


P43, i 


(A3) 


The resulting compatibility equation, which applies for tl 


nal) =P lh) + Poe 
""\4A, .E iGL/ j+l, i Vn jit (4 


Fic. A2. Details of panel. 


Subtracting Eq. (A3) from Eq. (A2) and noting that 


Gj. 1-1 — Qari = (P3i — Py, i-2)/Lia 
Gj,int ~~ Gein = (Pj, i+2 ite P;, ;) Lis 
and 
Hang —— Hyg = Xi 3 XY j-1, i Ww 
gives 


P ( 2L + =.) =P ( L ) 7 e 
ls) 0” eee Oe 41, i 


vg) + Paes (Ger) + Po (ier) 
P; i—2 Py i+2 = (A4) 
es si os tGL r tGL 


which, as indicated in Fig. A2, applies for all interior 
points except those next to the skewed rib. 

The compatibility relationship for the points adjacent 
to the skewed rib is derived by a procedure similar to 
that used to obtain Eq. (17). 
is somewhat different than the general expression be- 


The resulting expression 


cause of the following items: 
(a) Cap Energy Expression. 
cap energy term is taken as 


[(Po, 2 + P», 4) 2]? (Lo. 3 2A> 3E) 


A typical end bay 


Note that this expression differs from that used for an 
interior cap [see Eq. (7) ]. 

(b) Constraining Relations of Statics —The effects 
of the loads Ro... . Rss are included in writing the 
expressions of statics )> F = 0 and })M, = Oat rib sta- 
tions 2, 4,6, and 8. In addition, the relationship given 
by Eq. (A1) is used. 


1e interior points next to the skewed rib, is 


LV ma *,)4 
A jai, iii si A; 1,1 


w i, w 
F; ( _ ) F948 ( ) (A5) 
NGL” 44, aE) +o NGL 
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Eqs. (A4) and (A5) are then sufficient to determine the loads at any internal point in the panel in terms of the 


loads at adjacent points. 
and >> M = Oat each rib station. 
skewed boundary. 


Solution for Stringer Loads 


The external point loads are determined from the requirements of statics that }>F = () 
In addition, the requirement of Eq. (Al) must be satisfied for the points on the 
The foregoing equations are then sufficient to solve for loads throughout the panel. 


Evaluating the coefficients of Eqs. (A4) and (A5) term by term gives the following (note that these equations 


have been multiplied by G and that a value of G, /: 


0.385 has been used in the coefficient evaluation) : 


2 10 0.385 2 xX 10 LO X 0.4585 lO X 0.3885 
Py 4 ce - . -— re ) ’ Fy 4 ( a ) t P ssa, i ( - ) 
3.0 0.25 XK 10 2. 5.0 


Pi —_ 0.121 (P, 4 + P +1 


‘6. (° xX 10 X 0.385 _ 
j,t ‘ 
1 X 3.0 


Pe == 0.0294 Ugh +- Psat. 08 


These equations form the basis of the solution. 


The general approach for solution is a process of 
iteration. The details of the procedure are as follows: 
(1) Loads are assumed at all stations throughout 
the structure and are so selected that the requirements 
of statics are satisfied. Improved values of the loads 
are then obtained as follows: 
(2) 
(3) 
the free body obtained by cutting the structure just 
to the left of rib 4 and writing the equation that the 
sum of the moments about point 1, 4 (i.e., the point of 
intersection of stringer | with rib +4) is equal to zero. 
This equation is written as )0 1,4 = 0. 


P», 4is found from Eq. (A7). 
P,,2is found by considering the equilibrium of 


(4) Pi.4 is found from the consideration that 
> M3, 4 = 0. 
(5) P2,6is found from Eq. (A6). 
(6) Ps,61s found from Eq. (A7). 
(7) Ps.4 1s found from the consideration that 
> Mi.6 = 0. 
(8) Pi.6 is found from the consideration that 
~M,, 6 = U. 
(9) Ps:,sis found from Eq. (A6). 
(10) Ps, sis found from Eq. (A6). 
(11) Py, sis found from Eq. (A7). 
(12) Ps,6 is found from the consideration that 
> Ms = 0. 
(13) Pi,s is found from the consideration that 
> Mss = 0. 
(14) Ps, sis found from Eq. (A1), 
Ps.s = 4 — (Poo + Psa + Pas) 
(15) Loads at all other interior points are found by 


4) + 0.379(P;. «2 + Pj. i42) 
2 10 lO X 0.385 10 & 0.385 
+ IP) = as SM) + Pane (EO) 4 
0.25 X 10 t X 3.0 t X 3.0 
3S. 


10 X 0.385 10 
3 | Pe 
Py ( 3.0 ) il ce (—. x 10 


1) + O.118P 1, ¢ + 0.339P,, 1 + 0.367P), i441 


0) 


(A6) 


, ( 10 ‘ 10 
Pi 2 9- + f » 142 ae « 
XNO.25 X 10 O:20 


1O XK 0.585 
) Lp 


( 10 
1X 3.0 N\0.25 X 


0) 


(A7) 


Eq. (A6), and those at all exterior points are found by 
the requirements of statics. At the center of the panel 
span, the antisymmetry of the panel geometry should 
be kept in mind to determine the values for load in 
Eq. (A6) (i.e., Po, 1s must equal Py 14). 

This process is then repeated until the agreement 
between two successive passes is considered sufficiently 
close. In Fig. A3, the top numbers in each set represent 
the initial assumptions for loads. The second numbers 
are the estimates made after two passes, and the lower 
numbers are the final values obtained after four addi 
most 


tional passes. The last differences obtained in 


cases are of the order of 0.001 kips. 
ELECTRICAL ANALOGY SOLUTION 
Direct-Current Circuit 
The resistor and battery grid for the analogy is shown 
in Fig. A4. The magnitudes of the resistances in the 
circuit are Table 1. 
Since only the relative magnitudes of the resistances 


determined in accordance with 
need be maintained, the true values of resistance are 
multiplied by the factor 10.4/£ in order to obtain a 
convenient set of numerical values. 

In Fig. A4, the boundary currents, corresponding to 
the boundary shear flows in the structure, are known. 
The voltages @ .. . és, which are proportional to the 
Lagrangian multipliers XY. ... Xi in Eqs. (Al) and 
(A2), are initially unknown. 
determined such that the currents in the grid satisfy 


These voltages must be 


the requirements of statics in the actual structure. 
Solution for Currents and Voltages 
the 


The currents and voltages are determined by 


following iteration process: 
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Direct current analogy network for panel. 
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ONE AMPERE EQUALS ONE KIP. 
ARROWS INDICATE THE 


DIRECTION OF CURRENT 
IN THE CIRCUIT. 


0.900 
| UO) 


1.00 1.00 1.00 1.00 


Panel loads 


“oa 


(1) Currents that satisfy the 
statics in the structure are assumed throughout the net 
work. 

(2) Values of the voltages é . . . eis are obtained by 
applying Kirkhoff’s law [}>(/R + e) = 0] for each 
closed loop and solving for the corresponding value 
of e. Since the initial values for J are incorrect, the 
values obtained for e in each loop will vary across each 
bay. These values are averaged. 

(3) Using the average values of e, improved values 
of current are obtained by a procedure described in 
reference 4. 
summation around the loop is made, and a correcting 
current making this summation zero is added. These 
corrections are applied systematically throughout the 
grid. 

(4) After several passes (about four for this grid), 
the currents will approach values satisfying Kirkhoff’s 
law. However, the requirements of statics that > M = 
0 at each rib station will no longer be satisfied. To 


For each closed loop, the voltage drop 


correct this condition, the moment unbalance at each 
station is determined, and a set of correcting currents 
(In this 


problem, these currents were assumed to have a linear 


that reduce this unbalance to zero is added. 


distribution across the panel.) 


a | 199 
Onn ou rye —D 


requirements of 





1.00 } | 00 







1,00 1.00 1.00 1.00 


electrical analogy method 


(5) Steps 2, 3, and 4 are then repeated until the 


values converge. 


In Fig. A5, the top numbers in each set represent the 
initial assumptions for currents. The lower numbers 
are the values obtained after three complete correction 
cycles. It is seen that values are approaching the 
values obtained in the direct solution. In this example, 
it was found that solution of the equivalent electrical 
network by the procedure here described involved con- 


siderably more labor than the direct solution. 
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The Decay of Isotropic Temperature 
Fluctuations in an Isotropic Turbulence’ 


STANLEY CORRSIN?# 
The Johns Hopkins Uniwwersity 


SUMMARY 


The approach of von Karman and Howarth (to isotropic turbu 


decay of isotropic temperature 


has been applied to the 


lence 
fluctuations in a decaying isotropic turbulence, under the restric 


tion that the temperature differences are 


effect on the velocity field 
From the correlation equation it is found that, if the scalar func 
tions characteristic of the pertinent double and triple correlation 


ire assumed to descend to zero faster than r 2, the product of the 
r 
second moment of the double temperature correlation coefficient 


and the mean square temperature fluctuation is a constant dur 


ing decay his is analogous to Loitsiansky’ theorem” for 
the constancy of the product of turbulent energy and fourth mo 
turbulence de 


ment of the velocity correlation coefficient in the 


cay 
For extremely small Péclet 
to the three-dimensional heat conduction equation with 


Number the correlation equation 


reduces 


spherical symmetry Che corresponding equation for the linear 


viscous decrement of the velocity field is a quasi-five-dimensional 
equation 
growth of the temperature 


Péclet 


compared with the corresponding results for the 


Estimates of decay and scale 


fluctuations at both extremely small and extremely large 
Numbers are 
the most practical result is that the tem 
off at a 


Perhaps 


turbulence 


perature ““spottiness’” always seems to die lower rate 


than the velocity “spottines 


(1) INTRODUCTION 


ee THE BASIC PAPER of Taylor in 1935! intro 
ducing the concept of an isotropic turbulence and 
presenting the first contributions toward a statistical 
theory of this type of flow, a considerable amount of 
work has been done toward the development of a com 
plete theory. An exhaustive outline and reference list 
would be out of place here, but the introduction and use 
of correlation tensors by von Karman and Howarth, 
the demonstration of the Fourier transform relationship 
between correlation function and power spectrum by 
Taylor,* and the deduction of a part of the spectral (or 


von 


correlation) function by Kolmogoroff,4 Onsager, 


Weiszaicker,® and Heisenberg’ have all been major 


Much of the most recent work is contained in 
Karman*® and of 


steps. 


the publications of Lin and von 
Batchelor and Townsend.? 

The theoretical investigations of isotropic turbulence 
have proceeded along two partially independent lines 


Transfer and Fluid Mechanics 


30, 1950 
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Kovasznay and F 


too small to have any 


a Predictions of the time rate of change of the 
various characteristic statistical parameters during 
decay” under various simplifying restrictions. These 


parameters are usually the mean random kinetic energy, 


the ‘“‘scale’’ and the “‘microscale,’’ as defined by Taylor. ! 


(b) Derivations of the shape of the power spec 


trum (or correlation function) for various particular 


ranges of frequency or wave nuinber 
the work on (a 


The synthetic independence of and 


on (b), as well as the necessity for restrictions in (a 


and piecemeal treatment in (b), is a result of the ex 
treme difficulty of the problem, including the unfor 
tunate fact that in its statistical formulation the prob 
lem has more unknowns than equations. In spite of 
these hardships, attention to simple cases and the use 
of judicious postulates in the more complicated cases 
have been combined to yield a considerable body of 
theory on the statistical behavior of an isotropic, homo 
geneous, random velocity fluctuation field in an incom 
Furthermore, much of this 


pressible continuum. 


theory has been verified by experiment 
A closely 
ittention is that oi 


related problem that has apparently re 


ceived no the behavior of a similar 


calay field in a continuum. Obviously, the case of 


spottiness 
I 


simple (linear) diffusion of a random scalar 


is too special; an analysis of more general interest 


must include forced convection due to a concomitant 
The present paper gives part of an 
} 


nave proved 


field of turbulence. 
attempt to apply the various methods that 


inalysis 


useful in dealing with isotropic turbulence to an 
of the statistical behavior of an isotropic scalar fluctu 
ition field in an isotropic turbulence. Although the an 
alysis will be equally valid for concentration fluctu 


ations in the isotropic turbulent mixing between two 


gases of approximately equal molecular weights, it 


be set up as a temperature fluctuation theory it 


tempts at experimental verification can be done st 
easily on temperature fluctuation fields 

This paper is an introduction to the thermal problem 
analogous to (a)—1.e., attempts at prediction of the 


time rate of change of the significant statistical t 
parameters during simultaneous decay of an isotropi 
temperature fluctuation field and the isotropi 


lence in which it is assumed to occur 


are almost as much historical as logical, this aspect 
the problem has been handled in terms of correlatior 
functions rather than energy spectra. In fact, it does 


seem more convenient to discuss the behavior of a cd 
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caying fluctuation field from the point of view of the 
correlation equation and in terms of characteristic 
quantities (such as scale and microscale) which were 
first introduced in connection with correlations. How- 
ever, there is little doubt that the statistical behavior 
of a postulated stationary fluctuation field is more 
easily discussed in spectral terms. 

Some work has already been done on the shape of the 
temperature power spectrum in various wave number 
ranges. These results are published in the Journal of 
Applied Physics, March, 1951. 


(II) PRELIMINARY ‘‘KINEMATICS”’ 


For a temperature fluctuation field in an isotropic 
turbulence with zero mean velocity, the condition of 
isotropy is taken identical with that for turbulence 
alone: The average value of any function of the tem- 
perature fluctuations and their derivatives, defined rela- 
tive ‘to a particular set of axes, is unchanged by rota- 
tion or reflection of the axes. Such a definition ob- 
viously implies homogeneity of the field. 

The averaging is over space, although ensemble aver- 
ages over a large number of “‘universes’’ may be more 
satisfactory from a rigorous mathematical point of 
view. 

The analysis is restricted to temperature fluctuations 
sufficiently small that the density is effectively constant. 
Furthermore, the Mach Number is taken as approach- 
ing zero so that density changes due to pressure fluctu- 
ation are negligible, as is heat generated by viscous 


dissipation. 


Double Temperature Correlation 


In an infinite three-dimensional isotropic velocity and 
temperature fluctuation field, consider the correlation 
between instantaneous temperature fluctuations, 3 and 
dv’, at points P and P’, respectively. Because of isot- 
ropy, it is obvious that at any time this must be simply 
a scalar function of 7, the distance between P and P’. 


Hence, we introduce 


m(r, t) = dd’ 3, (¢ = time) (1) 


This is in contrast with the general velocity double cor- 
relation function, which is a second-rank tensor, though 
depending ultimately upon just a single scalar func- 
tion.” 

Since it depends only on the magnitude of PP’, m 
must be an even function of 7. Hence, when we ex- 
pand in series about 7 = 0, the coefficients of the odd 
powers must be zero. Also, from Eq. (1), it is obvious 
that m(0, t) = 1. Hence, for small 7, 


1 + m"(0, t) (r?/2!) + 
m” (0, t) (74/41) +... (2) 


mir, t) = 


The prime here denotes differentiation with respect to r. 
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Temperature Derivative Correlations 


An evaluation of correlations of type 
(Od /Ox;) (Od/Ox;) 


is of significance in a consideration of the conduc. 
tive ‘‘dissipation’’ of the temperature fluctuation field, 
We follow the approach used by von Karman and 
Howarth for turbulence: 

Let x; be the vector coordinate of point P and x,’ the 


vector coordinate of P’. Consider 


(0/Ox;) (30’) = 37(0m/Ox;) 
But also, | 
(0/Ox;) (88’) = 3'(08/O0x;) 
since J’ is independent of the choice of P. 
Equating these two expressions and differentiating 
with respect to x,’, we find 


(= \(e:) 92 O’m ‘ 
= —/ (3 
Ox; / \Ox,;’ OF OE 


where &, is the vector distance between P and P’, 


é = X;y . 


fi — Xe; eee = 7° 


Next, 7 is reduced to zero, so that 3’ ~ J. Then 
Eq. (3) can be written as 


( Ov \(2° ) m 
Ox; / \Ox; 


oe 1om E,£,0°m &,=, Om 
— 2° lim { 6, 
r—0) ro 7* Or y> OF 
Using Eq. (2), 
(Od /Ox;) (O8/Ox;) = —6,9?m"(0, £) (4 


6,, is the “Kronecker delta’’—i.e., equals 1 when i = 


j; equals 0 when? # j. In particular, when the index 
in Eq. (4) is repeated, 
(O8 /Ox;) (08 /Ox;) = —3d2m"(0, 6) (5 


Scale and Microscale 


It is convenient to follow Taylor' in introducing two 
characteristic quantities with the dimensions of length 


Lo = / m dr (6 
0 


The Thermal Microscale: 


9 
As” = 


The Thermal Scale: 


—2/m"(0, t) (7 


Clearly, \y is the r-intercept of the vertex-osculating 
parabola in m(r). 


With expression (7), Eq. (5) can be rewritten as 


(08 /dx,) (08 /Ox,) = 6 B2/ry? (8) 
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DECAY OF ISCGTROFPIC 


Double Temperature-Velocity Correlation 


yon Karman and Howarth? have shown that in iso- 
tropic turbulence the pressure fluctuation at a point is 
uncorrelated with the velocity fluctuation components 
at any point in the field. Since this proof is equally 
valid for the combination of any scalar fluctuation and 
the velocity, it follows that 


du; = 0 (9) 
u;' is the velocity fluctuation at P’. 


Triple Correlations 

The only type of triple correlation occurring explicitly 
in the present analysis is that between two temperature 
fluctuations and one velocity fluctuation, ui’, which 
we characterize by the dimensionless vector 


S, = u,00'/ Vu? #2 (10) 


The corresponding correlation in the turbulence is a 
tensor of third order. In order to deduce an expression 
for S; in terms of a scalar function of 7 and the vector 
distance £; between points P and P’, consider first 
the possible components of S; when one of the three 
orthogonal axes is set along the line connecting P and 
P’. By the required invariance to rotation and re- 
flection of the coordinate axes, it is obvious that the 
only nonzero component in this case is the one shown in 
Fig. 1. 

This component is a scalar function of 7, p(7), say. 
By inspection we see that it is an odd function of 7, 
since location of P’ on the —.x; axis would leave 1 point- 
ing away from P’ instead of toward it. By similar rea- 
soning it also follows that u,;/d0' = —u,dd’. 


Si(p, &) follows from inspection of the more general 
orientation, with P’ anywhere in the positive coordi- 
nate octant, at vector coordinate & = (£1, &, &3). 

From Fig. 2, by taking the projection of 1, along PP’ 


we see that 


ud’ (3? Vu? = cos B cos 46- p(7) 


Butcos 8 = £,/rcosé6. Therefore, 


u,00'/30? Vu? = pl(r) (&1/7) 


Similar examination of uw.’ and u;80’ leads to the gen- 
eral expression, for any time, 


S; = p(r) (&/7) (11) 


Since p(7, t) is an odd function of 7, its series expansion 
abour 7 = 0 contains only odd powers. Furthermore, 
the first derivative at 7 = 0 is zero, since, if we set the 


¥, axis along 7, 


re) 
[32 V u?]-! lim | 5 a0 | 


y—~e x1 


6'(0; i) = 


But, asx,’ x), 0’ ~> v8. Thus, 
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p’(0, t) ~ (oO Ox}) (1 ,3") = (0) 


due to homogeneity and to the zero correlation between 
u, and the scalar #*. 


Hence, for small 7, 
/ 
t * 9) 
+PGH.T--- 12) 


sve 


p(r, t) = p’’’(0, t) - 
3! 


(III) THE CORRELATION EQUATION 


Derivation 

From the incompressible heat-transfer equation, it is 
a simple matter to deduce an equation relating the 
double correlation function m(7, t) and scalar p(r, ¢) 
characterizing the triple correlation vector. 

We postulate a decaying ‘‘box’’ turbulence with no 
overall heat transfer. The heat-transfer equation for 
point P is 

Ob Ov a 0°08 = 
or Ox, pe, Ox Ox, | 
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where a = thermal conductivity and c, = specific heat 


at constant pressure. For point P’, 


Ov’ , OW . a O78’ 
“tt; 


= a (14) 
ot On; pC, OX;'OX; 


Multiplying Eq. (13) by 8’, Eq. (14) by #, averaging, 


and adding the two resulting equations, we get 





O ) 4 coe , 02( 32’) 
(3d) — (4, OF) — = (4 00") = 2y 
ol &; O&; * 0£,0&; 
where 
- . ad == . / - 
ef = = PCy; ry fn Xv, ~ Xx, 
Recalling that u,’33’ = —u,00', we combine the second 
and third terms and introduce the correlation coeffi- 
cients 
Oo OS, a: O’m 15 
(32m) — 2 8Vy2— = 27 077 (15) 
f ol & O0&,0§; 


With the introduction of Eq. (11), a little manipula- 
tion produces the final form of the correlation equation 


Op y) 
~ 29 Via (22 + 2°) = 
Or ? 


= (— =m) 7 
2yi9? (16) 
or? ae Orv 


eee oe ee - 


ete: 


(237712) 


oo 


This is similar in form to the well-known von Karman- 
Howarth equation for the velocity correlation func- 
The principal difference is the factor of 2 in the 
The 
term in Eq. (16) has the form of the 


tions. 
lower order terms in the two large brackets. 
“dissipation” 
spherically symmetrical Laplacian in three dimensions, 


We choose 6 = 2. Then, if P(r), m(r), and Om/Or de- 
crease to zero faster than 7~? as 7 > o, the right-hand 


side is zero. Hence, during decay, 


ye / mr? dr = const. = N 


aw @ 


From an examination of the power spectrum of the 
temperature fluctuations at extremely low wave num- 
bers, the suitable rate of decrease of the correlations at 


large 7 can be demonstrated. 


AEBRONATTTICAL 
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while in the velocity correlation equation, as pointed 
out by Loitsiansky,'"’ the corresponding term has the 
form of a spherically symmetrical Laplacian in five di. 
mensions. 

The decay equation is obtained from Eq. (16) by 
taking the limit form as 7 > 0: 


dg?/dt = Bbyd2m" (0, t) (17 


or, in terms of thermal microscale, 


dg?/dt = —12y(82/X,") (18 


which is just like Taylor’s turbulence decay equation 


dy?/dt = —10v(u2/r?) (19 


where v = kinematic viscosity and \ velocity micro- 
scale corresponding to the vertex-osculating parabola 


for the von Karman-Howarth g(7, ¢) correlation function 


An Invariant During Decay 


Loitsiansky has shown" that if the double and triple 


velocity correlations in turbulence (and_ their first 


derivatives) decrease more rapidly than 7~‘ as 7 > 
then there is a simple invariant during the decay of 


turbulence—namely, 


wf r*f dr = M = const. (20 
we 


where f is the scalar double velocity correlation for which 
the velocity components are taken along the line con- 
necting P and P’. 

An analogous invariant can be deduced for tem 


perature fluctuations. We follow the procedure of 


Loitsiansky, multiplying Eq. (16) by 7 and integrating 


over 7 from zero to infinity. There results 


7" ; - 7s c b om , 
(wf mr’ ar) = 2? V/42| pr’ + (4 — 2b)8?~ wf pr’ * dr + 20° °| } 
dt 0 0 0 or 0 


*? Om 
2¥(2 — byw? yo dy 
0 Or 


(IV) Decay OF TEMPERATURE FLUCTUATIONS 


Case of Extremely Small Péclet Number* 


A Péclet number characterizing the combined heat | 


and velocity system may be chosen as 
2 —_ > A 
F oo u2r» Y 


When this becomes sufficiently small, the triple cor- 
relation term in the correlation equation may be neg- 


* It has been called to my attention by Dr. C. C. Lin that 
this degenerate case is somewhat like the one worked out by 
Reissner (Proc. 5th Int. Cong. Appl. Mech., 1938) for (erroneous 
application to the velocity field. 
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DECAY OF ISOTROPIC 
lected so that 


9 9 


- (92m) + - 
Oo y 4 


vy yr O7 


(32m) = 2y - (3m) (22) 
where #2 is put inside the 7-derivatives for simplicity of 
form. We are now dealing with a simple (linear) 
random heat-conduction problem in three dimensions, 
and it is only slightly surprising that Eq. (22) should 
turn out to be the three-dimensional heat-conduction 
equation with spherical symmetry. The double tem- 
perature correlation function is, after all, spherically 
symmetric in the physical space—unlike the double 
Recalling that the linear 


(22) 


velocity correlation tensor. 
velocity decay equation corresponding to Eq. 
has the form of a five-dimensional spherically symmet- 
ric heat conduction equation, we observe that in some 
way the extra “‘degrees of freedom’’ in the dynamical 
problem have the same effect as would extra dimensions 
in the physical space. 
In any case, the solution to Eq. (22) which corre 
sponds to that used in the turbulence decay is 
v2m(7, t) const. pos 3 =v 8y7)] (23) 
(yr) 


where 7 f f;, is the time variable. /; is the time 
between the apparent origin of the thermal fluctuations 
and the apparent origin of the velocity fluctuations, 
arbitrarily chosen at ¢ = 0. 

The turbulence decay for negligible inertia effects 
small Reynolds Number) was given completely by 


Loitsiansky aa 
Svt)| 


M_ exp [-— (7 


INa/2r (vt) 


u2f(7, t) (24) 


The integration constant in Eq. (23) can be evaluated 
from the invariance relation (21) so that 


N exp[— (7?/8yr) | 
v?m(r7, t) = ( ) : : (25 
14/27 (yr) °° 


The time variations of 32, Ay, and Ly are easily ob- 


tained from Eq. (25) as 


v2? = (N/4V 2m) (y7) (26) 
Ay” SYT (27) 
Ls? = 2ryr (28) 


These are to be compared with the turbulence decay 


results: 
u2 = (M/48vV/2r) (vt) (29) 
A? = 4p (30) 
L? = (w/2)vt (31) 


where both \ and Z pertain to the von Karman-Howarth 
g(r, t) double correlation, with velocity components at 
P and P’ taken perpendicular to the line connecting the 
points. 

The principal conclusion of a comparison between 
these two sets of results is that, when both Reynolds 
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Number and Péclet Number are extremely small, the 
temperature “‘spottiness’’ dies out at a slower rate than 


the velocity “‘spottiness.”’ 
For times sufficiently large so that t; < t, we have 
7 = t, Then the thermal and velocity scales are re 


lated by 


Ay?/A? = 2/0 (32) 
L,?/L? = 4/¢ (33) 


where o uc,/a, the Prandtl Number. The numer- 
ical constants in Eqs. (32) and (33) are not completely 
indicative of the relative magnitudes of the phenom- 
ena; the velocity system also has scale L, = 2L and 
microscale A, = /2A which are characteristic of the 
von Karman-Howarth f(7, t) correlation. Hence, in- 


? 


stead of Eqs. (32) and (33), we could write 


for these extremely ‘‘tired”’ fluctuations. 
To conclude this section, we observe that, at small 
Péclet Number, 


Page (36 
while at small Reynolds Number 
R,~ tit” 37 


The Reynolds Number definition is 


R, V u?r/v 


Case of Extremely Large Péclet Number with Self- 
Preserving Correlations 


Next, we follow the lead of von Karman? in looking 
for information on the time history of the characteristic 
parameters when the triple correlation terms dominate, 
under the restricted 


postulate of ‘“‘self-preserving”’ 


correlation functions. By self-preservation we mean 


here that 


m= mQ); p= P(2 38 
where Q= 7/L,. 

In this case we use Ly as primary length rather than 
Ay, Chiefly because conduction plays a relatively in- 
direct role as contrasted with the previous case. Hence, 
we introduce a characteristic Péclet Number, 


Pi = V yrho/1 

analogous to the Reynolds Number 

R, = Vy2L/v 
Introducing assumption (38) into the correlation 
equation (16), assuming that P, is sufficiently large to 
permit neglecting the right-hand side and replacing 


dj?/dt by its exact expression from Eq. (18), we end up 
with the approximate equation 
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19 Le (9) + dLy " dm oe dp 4 *) ‘ 
2 m Q 2./-3 2 = 
Lae a dat 7V wag Q 


(39) 
Our restriction of self-preserving m and p can be 
satisfied only if the three coefficients, which are func- 
tions only of ¢, are proportional. 
This condition gives two equations for the charac- 
teristic variables: 
(40) 
(41) 


ur," yLy = c 
dL,/dt = DV/q2 


where C and D are constants. These two equations 
correspond to the pair deduced by von Karman for 
the turbulence: 


u2r*?/vL = A (42) 
dL/dt = BV 42 (43) 
where A and B are constants. 
Dividing Eq. (41) by Eq. (43), we get simply 
dL,/dL = D/B (44) 


and with the boundary conditions that at ¢ = ft, L = 
Lo, and Ls = Lo, 


Ly Dly (L . 
= ] + — | (45) 
Le B | ee Lo 
with von Karman’s result for L(¢)—1.e., 
s Lo = (t = sil ae (46)* 


we obtain 


* This result and some of the following ones from reference 2 


were misprinted there. 


oO ~ capt = 


| 12 ) if 
CD\i5 + AB 


This compares with 


= ( A yey" (6 +4) 
- - 5 + AB to to 


kT . c 1B/(5 1B) 
Therefore we re-examine the problem for (t/fy)*”” ° * 


a i. Then, if Ly,/Zo and D/B are both of order 
unity, 
L/L = D/B = const. (55) 
Ae*/A* = (CD/AB) (1/c) (56) 
and 
ce Dy? 4 (t care CD)[12/(5 + AB) ] (57) 


where #,2 is the value of 3? at some specified time fo, and 
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Ly D Ly 


, AB/(5 + AB) 
wai l ( 7) 
Ln Pe. |(;) | a 


Thus, although Ly is not simply proportional to L, 
Eq. (44) tells us that it has the same rate of growth, 
The ratio of the two scales is 


Ls me D 4 (> x~ (54+ AB) 
LB Lo B/\to (48) 


This approaches a constant only for extremely large ¢, 
Division of Eq. (40) by Eq. (42) gives the micro- 
scale ratio in terms of the scale ratio: 


Xg?/A* = (C/A) (1/e) (Ly /L) (49) 
von Karman found 
A? = o(5+ AB)t (50) 
so that 
A»? . (5+ AB |p +(+ *)(.) a ae 
ag? : AB) = 
PS“ "se '°\i, ore 
(51) 


In passing, we note the analogous relations for the 
ratios of microscales to scales, obtained directly from 
Eqs. (40) and (42): 


‘ike + i/P~ 1/77, (52) 
and 
(53) 


A\/L ~ 1/R, ~ 1/VR;, 


To find 3*(¢), we substitute \, into the decay equation 
[Eq. (18) ]. 
plicated for simple physical inspection : 


This yields an expression that is too com- 


dt/t 


Di (2 > (*) ied ds 
B Lo B lo 


é 12 B Ly, 
K = exp | = ; (1 - )| 
CD D Lo 


which, to be consistent with the assumption of this 
special case, is of order unity. 

Eq. (57) may be compared with von Karman’s re- 
sult for the large Reynolds Number (self-preserving) 
turbulence at any /: 


9 » — (5 + AB on 
u?/uo2 = (t/to)~*” ss (58) 


The relative rates of decay of velocity and temperature 
fluctuations remain unspecified until estimates are 
made of AB and CD. 

The results in Eqs. (55), (56), and (57) are valid only 
with extremely large Péclet Number and (t/f)“""° * ia 
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> 1. This is a difficult pair of conditions to meet 
simultaneously. On the other hand, virtually identical 
relations hold at all values of (¢//)) if only it is possible 
to make Ly/L = D/B at time f. Presumably, in a 
wind tunnel, Ly and Ly, can be set independently if 
the turbulence is produced by a grid set in a uniform 
air stream and the temperature field is produced by a 
square-mesh lattice of heated wires whose Reynolds 
Number is too small to produce additional turbulence. 
For all (¢/¢o) in this particular case, Eqs. (55) and (56) 
hold, while Eq. (57) is replaced by the simpler form 
82/9,2 = (t/to) ~ AP/CPM/E + 4) (59) 
In these two special cases, the time behavior of the 
Péclet Number is identical with that of the Reynolds 
Number in the large Reynolds Number turbulence 
at any (¢/to): 


Pye GPO F® (60) 

This leads us to suspect that any reasonable estimate 
should give AB < 5. 

The final step in this phase of the analysis is the 

evaluation of CD from the invariant relation (21). 

With complete self-preservation, Eq. (21) implies that 


§2L,? = const. (61) 


whereas the Loitsiansky relation (20) implies 


nu? L® = const. (62) 

Kolmogoroff'! showed that, with Eq. (62), AB = 2, 
so that 

L/Lo = (t/to) (63) 

2 = 7pt (64) 

u?/Uy? = (t/to)~” (65) 


for any (t/to). 

For the temperature fluctuations, the application 
of Eq. (61) to the general, self-preserving, large Péclet 
Number system is too complicated for the crude aims of 
this analysis, so we again inspect only the special case 
with Ly,/Lo set equal to DB. Then Eq. (61) gives 
CD = 4, so that 


Lo/Lo, = (t/to)’ (66) 
Ay? = 14yt (67) 
32/9,2 = (t/to)~ ” (68) 


A comparison of these relations with the turbulence 
results above shows two estimates of particular in- 
terest: First, 

A,?2/k2 = 2/e (69) 


just as in the case of extremely small Reynolds Num- 
ber and Péclet Number when ¢ > ¢;. The second 


to 
w 


point is that here, too, the temperature fluctuations die 
out at a slower rate than the velocity fluctuations. 

it should be noted that the general problem postu- 
lated in this paper can be easily set up in the laboratory 
by using a heated grid to produce isotropic turbulence. 
This case would presumably fix L,; ~ L. Furthermore, 
as implied earlier, the experiments can be set up for 
widely varying L,/L if the temperature field is pro- 
duced by a fine-wire square-mesh lattice rather than 
by the turbulence grid. Of course, all pertinent sta- 
tistical quantities in both dynamical and thermal sys- 
tem can be measured with hot-wire anemometry equip- 


ment. 


REFERENCES 


' Taylor, G. I., Statistical Theory of Turbulence, Proc. Roy 
Soc. London (A) (No. 873), Vol. 151, pp. 421-478, September 2, 
1935. 

2 von Karman, Theodore, and Howarth, L., On the Statistical 
Theory of Isotropic Turbulence, Proc. Roy. Soc. London (A) (No. 
917), Vol. 164, pp. 192-215, January 21, 1938 

3 Taylor, G. I., The Spectrum of Turbulence, Proc. Roy. Soc. 
London (A) (No. 919), Vol. 164, pp. 476-490, February 18, 1938. 

‘ Kolmogoroff, A. N., The Local Structure of Turbulence in In- 
compressible Viscous Fluid for Very Large Reynolds Numbers, 
Comp. Rend. (Doklady) de l’'Acad. Sci. URSS, Vol. 30, No. 4, 
pp. 301-305, 1941. Also, Dissipation of Energy in the Locally 
Isotropic Turbulence, Comp. Rend. (Doklady) de 1l’Acad. Sci. 
URSS, Vol. 32, No. 1, pp. 16-18, 1941. 

5 Onsager, L., The Distribution of Energy in Turbulence, 
Physical Review, Vol. 68, Nos. 11 and 12, p. 286, December 1 
and 15, 1945 

6 yon Weiszacker, C. F., Das Spektrum der Turbulenz bei Grossen 
R-ynoldsschen Zahlen, Zeits. Phys., Bd. 124 (Heft 7/12), pp 
614-627, September &, 1948 

7 Heisenberg, W., Zur Statistischen Theorie der Turbulenz, Zeits. 
Phys., Bd. 124 (Heft 7/12), pp. 628-657, September 8, 1948. 

® Lin, C. C., On the Law of Decay and the Spectrum of Isotropic 
Turbulence, Proc. 7th Int. Cong. for Appl. Mech., 1948. Also, 
Lin, C. C., and von Karman, Theodore, On the Concept of Simi- 
larity in the Theory of Isotropic Turbulence, Rev. Mod. Phys., 
Vol. 21, No. 3, pp. 516-519, July, 1949; Lin, C. C., Remarks on 
the Spectrum of Turbulence, Proc. 1st Symp. in Appl. Math 
(1947), published by Am. Math. Soc., 1949; etc. 

9 Batchelor, G. K., and Townsend, A. A., Decay of Isotropic 
Turbulence in the Initial Period, Proc. Roy. Soc. London (A) (No 
1035), Vol. 193, pp. 539-558, July 21, 1948. Also, Decay of 
Turbulence in the Final Period, Proc. Roy. Soc. London (A) (No 
1039), Vol. 194, pp. 527-548, November 9, 1948; The Nature 
of Turbulent Motion at Large Wave-Numbers, Proc. Roy. Soc 
London (A) (No. 1057), Vol. 199, pp. 238-255, October 25, 1949; 
Batchelor, G. K., The Role of Big Eddies in Homogeneous Turbu- 
lence, Proc. Roy. Soc. London (A) (No. 1043), Vol. 195, pp. 513 
532, February 3, 1949; etc. 

10 Loitsiansky, L. G., Some Basic Laws of Isotropic Turbulent 
Flow, Central Aero-Hydrodynamic of Inst., Moscow, Rept. No 
440, 1939. (Translated as N.A.C.A. T.M. No. 1079.) 

1! Kolmogoroff, A. N., On Degeneration of Isotropic Turbulence 
in an Incompressible Viscous Liquid, Comp. Rend. (Doklady) de 
l’Acad. Sci. URSS, Vol. 31, No. 6, pp. 538-540, 1941. 











PF 25028 F ES Fes 


= 2se 


Note on Bending of Thick Sandwich Plates’ 


GEORGE GERARD? 
New York Unwersity 


SUMMARY 
A system of equations describing the bending of thick isotropic 
sandwich plates is developed. The thick-plate system is re- 
quired when a characteristic length such as a buckle wave length 
is of the order of magnitude of the sandwich thickness. When 
the characteristic length is large, the thick-plate system reduces 
to the ordinary thin-plate system. 


INTRODUCTION 


— THEORIES for homogeneous plates are of 
three types, which essentially depend upon the 
thickness of the plate, /, relative to either the deflec- 
tion, 6, or a characteristic length, /: 

(a) Linear bending theory for thin plates is used 
when the deflections and thickness are small—t.e., 6/h 
<landh/l<1. 

(b) Nonlinear or large deflection theory for bending 
of thin plates is required when the thickness is small, 
h/l < 1, and the deflections are of the order of 
magnitude of the thickness. 

(c) Bending theory for moderately thick plates is 
based on the assumption that the deflections are small, 
6/h < 1, but that the thickness relative to the char- 
acteristic length is not small; hence, this theory in- 
cludes the deflections due to shear. 

A theory for sandwich plates corresponding to the 
conditions specified for Case (a) for homogeneous plates 
has been given by Libove and Batdorf,' among others, 
and also by an order of magnitude analysis’ that de- 
pends upon the conditions 6/h < 1 and h/l < 1 
for the sandwich plate. Since core materials are 
usually weak in shear, these theories necessarily account 
for the deflections due to shear and consequently em- 
brace the theory of Case (c) also. Corresponding to 
Case (b) for the homogeneous plate, a large deflection 
theory for sandwich plates has been given by Reissner.* 

The need for a third type of theory for sandwich 
plates arose in connection with an investigation of the 
instability behavior of thick sandwich plates for which 
wrinkling is the usual mode of instability. In this case, 
the deflection relative to either the face or sandwich 
thicknesses can be assumed to be small, but, since 
wrinkling is a short wave-length phenomenon, the 
thickness of the sandwich may be of the same order of 
magnitude as the characteristic length. 

Received October 5, 1950. 

* Sponsored by the Office of Naval Research under Contract 
No. N6-onr-279, Task Order V. 

+ Assistant Professor of Aeronautical Engineering, Guggenheim 
School of Aeronautics, College of Engineering. 


In this note, a general system of equations for the 
bending theory of ‘‘thick’’ sandwich plates is developed 
using a method introduced by Reissner* for thin sand- 
wich plates and incorporating a more complete descrip. 
tion of the core behavior such as used by Williams, 
Leggett, and Hopkins.‘ It is shown that the thick. 
plate system reduces to the thin-plate system when the 
characteristic length becomes much larger than the 
sandwich thickness—i.e., /1]< 1. 

In this development, it is assumed that the core and 
face materials are isotropic and that, although the core 
is thick, the faces are still thin relative to the char- 
acteristic length so as to behave as thin homogeneous 
plates. As stated previously, the need for this theory 
arose in connection with the wrinkling problem, but it 
may have further applications in various bending 
problems for sandwich plates. 


SYMBOLS 


D; = bending rigidity of the face about its middle sur- 
face, D; = E,t?/12(1 — v?,) 
e = average strain of sandwich in z-direction, 
e = (w, —wWe)/h 
E = modulus of elasticity 
G = shear modulus, G = E/2(1 + v) 
h = core thickness 
h = centroidal distance between faces, h = h + t 
l = length of plate 
M = moment resultant 
N = force resultant 


= transverse pressure 


@) = shear resultant 

t = face thickness 

u,v, Ww = displacements 

Ww = average displacement of sandwich, w = (w, + w»)/2 
x, ¥, Z = coordinate axes 

o = normal stress 

T = shear stress 

v = Poisson’s ratio 

r = m/l 

0 = \h/2 


Quantities with a bar denote core parameters 


Subscripts 
i= = upper or lower face (Fig. 1) 
ri = face parameters 


THICK-PLATE SYSTEM 


For simplicity, the sandwich plate is assumed to act 
under a condition of plane strain. Thus, the stress- 
strain relations for the faces are given by 

Ey Ou; 
: 92 (1) 


Oy = e=1,2 


(1 — vf) Ow’ 
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where the subscript 7 refers to the upper or lower face 
as shown in Fig. 1. In addition, the usual assumptions 
concerning the z-component strains that are introduced 
in deriving linear plate bending theory apply to the 
faces. 

The equilibrium equations for the faces follow from 
plate bending theory with the introduction of addi- 


tional terms where necessary as shown in Fig. 1. For 
the upper face, 
(ONy/Ox) — F221 = 0 (3) 
D,(0'w,/Ox') — (pi/2) + 6a = O (4) 
and for the lower face, 
(ONz2/0xX) + Tr2 = O (5) 
D(O4wW2/Ox*) + (p2/2) — G2 = O (6) 


In Eqs. (4) and (6), it is assumed that the core is at- 
tached to the faces at the middle surface of the latter. 
With this assumption, an additional moment con- 
tributed by the shear force is neglected, which is a 
reasonable assumption for thin faces. 

By the addition and subtraction of Eqs. (3) and (5) 
with Eqs. (4) and (6), the following combined equilib- 


rium equations are obtained : 


(ON,,/OX) — (F221 7-2) = 0 (7) 

(OM,,/O0x) — [A(Fr2n + 7r2)/2] = 0 (8) 
2D,(0'w/Ox*) — [(pr _ p2)/2] + (é4 — G2) = 0 (9) 
hD(d%e/Ox*) — [(pi + po)/2] + (@a + G2) = O (10) 


SANDWICH PLATES 


where 


Ney = Nu + Na 
M,, = (Na — N,2)h/2 


From Eq. (1), by addition and subtraction of the 
equation for each face, the following distortion equa- 
tions are obtained: 


Ni = [Ept/(1 — v/*)] [O(a + au) /dx] (11) 


M,, = [Eyth/2(1 — vf)] [O(a — u2)/Ox] (12) 

Although the faces are taken to act under a condi- 
tion of plain strain, a further assumption that Poisson's 
ratio for the core is zero is in reasonable agreement with 
many core materials and serves to simplify the equa- 
tions that describe the behavior of the core. With 
this assumption, the stress-strain relations for the core 


reduce to 


¢, = E(0u/Ox) (13) 
&, = E(Ow/dz) (14) 
712 = G[(Ou/dz) + (Ow/dx) | (15) 


The solution for the core stresses reduces to finding 
a suitable function ¢ which is a solution to the com- 


patability equation 


(0'¢/O0x*) + 2(04¢/0x70z7) + (O4~/02*) = O (16) 
The stress function can be taken as? 
gy = sin Ax(C; cosh Az + Cy sinh Az + 
C3z cosh Az + Cyz sinh Az) (17) 
The core stresses are 
a, = 0*¢/02" (18) 
&, = 0°9/O0x" (19) 
Tr2 = —O*e/O0x02 (20) 


The core displacements are obtained by integration 
of Eqs. (13) to (15) and use of Eqs. (18) to (20). 


—(1/E) cos \x[C,A cosh Az + Cod sinh Az + 


“= 
C;(2 sinh Az + Az cosh Az) + 

C,(2 cosh \z + Az sinh Az)] (21) 
w = —(1/E) sin Ax[C)A sinh Az + CA cosh Az + 


C3(Az sinh Az — cosh Az) + 
C,(Az cosh Az — sinh Az)] (22) 


The independent equations for the faces and the core 
can be combined to obtain the characteristics of the 
sandwich element by satisfying the continuity require- 
ments at points where the faces are assumed to be at- 
tached to the core—i.e., the middle surface of the face 
(see Fig. 1). These requirements at z = +//2 are: 
(a) the displacements u and w must be continuous, 
(b) the derivatives of « must be continuous, and (c) 


the stresses ¢, and 7,, must be continuous. 
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From continuity requirement (a), it is found that 


U, + ue = —(2/E) cos rx X 
[C\A cosh 6 + C,(2 cosh 6 + @ sinh @)] (23) 


Uy — Uz = —(2/E) cos dx X 
[C2 sinh 6 + C;(2 sinh 6 + @cosh @)] (24) 


From continuity requirement (b) and by substitution of 
Eqs. (23) and (25) into Eqs. (11) and (12), 


2k jt sin Ax 


N., = 
‘(1 = »PE 
[C:A? cosh @ + C,d (2 cosh 6 + @ sinh @)] (25) 
Eyth sin Xx 
M,y = 


(1 — VP)E 
[C.d? sinh 6 + C3 (2 sinh 6 + @cosh @)] (26) 


Similarly, from the definitions of the average dis- 
placements and Eq. (22), it is found that 


w= —(1/E) sindvAx X 
[CA cosh 6 + C;(@ sinh @ — cosh @)] (27) 


e = —(2/Eh) sin \x X 
[C\A sinh @ + C,(@ cosh 6 — sinh @)] (28) 


For completeness, the following additional terms 
that appear in the equilibrium equations are obtained 
by observing continuity requirement (c). 


Fra + T7122 = —2A cos Ax X 
[C:\ cosh 6 + C;(cosh 6 + @ sinh 6)] (29) 


Tr2l — Tr22 = —2) cos Ax b 4 
[C\A sinh 6 + C,(sinh @ + @ cosh @)] (30) 


G1 + G62 = —2d sin Ax(C\A cosh 6 + C,@ sinh 6) (31) 
G1 — G2 = —2d sin Ax(CA sinh 6 + C36 cosh 6) (32) 


Equilibrium Eqs. (7) to (10) and Eqs. (25) to (32) 
constitute a complete system of equations for the sand- 
wich plate under a condition of plane strain. For the 
solution of bending problems, it is necessary to deter- 
mine the four constants that appear in Eqs. (25) to (32) 
from boundary conditions at z = +(h/2 + 2). 


In the solution of buckling problems, the lateral 
loads p; and 2 vanish and the z-component membrane 
forces that are no longer negligible are introduced. In 
such problems, a method introduced by Williams, 
Leggett, and Hopkins‘ may be used in which the condi- 
tion of nonvanishing of the constants is used to deter- 
mine the buckling load. 


THIN-PLATE SYSTEM 


One should expect the analysis presented herein for 
thick plates to yield the essential results of the thin- 
plate analyses'~* by letting 


6 = dh/2 = mah/21 0 (33) 


AERONAU 
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In the thin-plate system, the face parallel stresses jn 
the core are neglected so that, from the equilibrium 


equation, 


(OG, Ox) + (OFxy Oy) + (O7;- Z) = 0 (34) 
the result obtained is 
07,,/02 = 0, Trz = f(x) (35) 


From use of Eq. (35), further results that characterize 
the thin-plate systems! ~* are obtained. 


Q, = hi,, (36) 
(Tra Tr2) = (G21 — G22) = 0 (37) 
(1/2)(é1 + G2) = Ee (38) 


To obtain the equivalent of Eq. (36) from the present 


analysis, let 
"h 2 
Q, = Tre OZ (39) 
e h/2 


By substituting Eq. (20) into Eq. (39) and integrating, 
we obtain 


QO, = —2 cos A x(CA sinh 86 + C36 cosh 8) (40) 
When @ is small, then sinh @=@andcosh@=1. Thus, 
lim QO, = —h(CoA + C3)A cos Ax (41) 
6—0 
From Eq. (29), however, 
lim Ce + T1292) = —2(Cor + C3)X cos AX (42) 


6—0 
By combining Eqs. (41) and (42), 


lim Q, = A(722 + F22)/2 
6—0 
which is the essential result embodied in Eq. (35). 
From Eqs. (30) and (32) we obtain results correspond- 
ing to Eq. (37). 


lim (7,21 bas Tr22) = —2Xr COs Ax( Cid + 2C,0) > 0 
6—~0 

lim (6: — G2) = —2dA sin Ax(C.A\8 + C36) — O 
a@—0 


Results corresponding to Eq. (38) are obtained as 
follows. From Eq. (28) 


lim e = —(2C,\0 Eh) sin Ax = —(C,\?/E) sin Ax (43) 


6—0 
and from Eq. (31), 


lim (G1 + 62) = —2d sin Ax(C\A + C6?) = 
aioe —2C,A? sin Ax (44) 


Therefore, from Eqs. (43) and (44), in the limit as 
Gq, 


(1/2)(@1 + G2) = Ee 


which corresponds to Eq. (38). 


(Continued on page 432) 
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RIEF REPORTS of investigations in the aeronautical sciences 

and discussions of papers published in the JOURNAL 

are presented in this special department. The publication will be 

completed approximately 8 to 10 weeks after receipt of the material 

The Editorial Committee does not hold itself responsible for the 

opinions expressed by the correspondents. Contributions should 
not exceed 800 words in length. 


On the Nonexistence of Subsonic Conical 
Flows 


Theodore R. Goodman 
Cornell Aeronautical Laboratory, Inc., Buffalo, N.Y. 
March 15, 1951 ‘ 


— PAPER! ON SUBSONIC CONICAL FLOWS, which appeared in 
the February, 1951, issue of the JoURNAL, seems to require 
some clarification. 

When the theory of supersonic finite wings was young, the 
early investigators tried to extend the known methods of sub- 
sonic wing theory to supersonic wing problems (e.g., Schlichting). 
Now the picture seems to be reversed—i.e., a known method of 
supersonic wing theory is being applied to subsonic problems. 
Although subsonic conical solutions may be used to solve non- 
conical problems, still there exist no subsonic conical solutions 
that satisfy both the boundary condition on the wing and at in- 
finity. This paradox needs further explanation. 

For ease of visualization, picture a delta wing of infinite extent 
at zero yaw flying at a subsonic speed. If a conical solution exists, 


the perturbation velocities are constant along rays passing 
through the apex of the wing, and, since the perturbation veloci 
ties are zero far ahead of the wing (boundary condition at in- 
finity ), it follows that they must be zero everywhere upstream of 
the plane that passes through the apex of the wing and is normal 
to the free-stream velocity vector. In other words, there is a 
jump discontinuity in the 


derivative of them at this plane or at some surface that extends 


perturbation velocities or some 


downstream of this plane. But if this is so, the velocities 
cannot be harmonic functions of the space variables, which 
they must be in order to satisfy the complete boundary- 


value problem. It necessarily follows, therefore, that no sub- 
sonic conical solution exists which satisfies the boundary condi- 
tion at infinity. For supersonic flow, on the other hand, the 
jump discontinuity is correct, and, in fact, the velocities are not 
harmonic functions but, instead, they satisfy the wave equation. 

Examination of the conical solutions obtained by Mr. Lampert 
discloses the fact that they do not satisfy the boundary condi- 
However, when two such solutions are combined 


as in the case of the infinite swept 


tion at infinity. 
to form a nonconical solution 
wing—the result does, indeed, satisfy the correct boundary con- 
dition at infinity. 
mathematical tool for solving nonconical flow problems, and it 
must not be construed that they have any significance in them- 


The conical flows, therefore, are merely a 


selves. The efficacy of using such a procedure instead of an 
alternative method (integrating doublets over the plan form, say ) 


must be weighed in any particular case. 


REFERENCE 


Seymour, Conical Flow Methods Applied to Uniformly Loaded 
18, No 


Lampert 
Wings in Subsonic Flow, Journal of the Aeronautical Sciences, Vol 
2 1951. 


<. p. 107, February 


On the Calculation of the Rate of 
Heat Transfer Through a Laminar 
Boundary Layer* 


Arthur N. Tifford and John Wolansky 

Aeronautical Engineering Department, The Ohio State University, 
Columbus, Ohio 

March 5, 1951 


D URING RECENT MONTHS, three general, relatively precise, 
approximate methods of analysis of the rate of heat trans 
fer through a laminar thermal boundary layer have been de 
scribed in the JoURNAL;'~* namely, Squire’s method,‘ Diene- 
mann’s extensiont of the Holstein-Bohlen-Pohlhausen method,? 
and a modification of Lighthill’s method.* The results obtained 
by means of the first two methods have been compared with the 
results of Goland’s more exact analysis in the case of a particular 
isothermal cylinder. It seems appropriate that the result of the 
application of the third method should be similarly compared. 
Fig. 1 gives such a comparison. 

Curves A and B of Fig. 1 represent applications of the modified 
Curve A has been determined by means of the 
In determining Curve 


Lighthill method. 
procedure recommended in reference 3. 
B, however, the procedure was altered to the extent of using an 
exact determination of the momentum thickness at each point. 
The two curves are significantly different only in the adverse 
pressure gradient region beyond the point of minimum pressure. 
The difference is explained by the rather poor approximation of 
the momentum thickness given by Falkner’s graph in that re- 
gion. Even so, it might be noted that the result of the more 
approximate analysis, Curve A, is in better agreement with Go- 
land’s reference cutve than are the results of his application of 
Squire’s analysis. Of course, the more exact Curve B agrees 
even more closely with the reference curve 

It is interesting to notice that Dienemann’s results and those 
of the present note are in essentially exact agreement up to the 
point of minimum pressure. Furthermore, Squire’s method of 
analysis yields heat-transfer coefficients in good agreement with 
these results in the region near the point of minimum pressure 
Goland’s so-called exact curve, on the other hand, yields heat- 
transfer coefficients some 5 per cent lower in this same region. 
Since the three independent approximative determinations should 
be particularly precise there and are in agreement, the precision 
of Goland’s exact curve appears to be suspect. The precision of 
the evaluation of the functions ot reference 6 upon which Goland’s 
values depend are, in turn, subject to some question. 


* Work supported by the Office of Air Research, U.S.A.F 
+ Seban has reported a similar analysis in reference 5 
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It does not seem out of order to report here that the reference 
exact analysis is identical with the first portion of a more com- 
plete general thermal boundary-layer solution obtained by the 
senior author of this note in some unpublished work in 1943. 
More specifically, it was shown that there is associated with a 
local free-stream velocity variation, such as that of Eq. (1), 


uy + aix + asx? + asx® + ayx7 +... Ct)” 
a thermal boundary-layer solution, Eq. (2), 
T- 7; 1 ’ a3” 
a — =] xd, + agx®d;s + x1 ass + P55 } + 
Ti = Tw uy qa 
i 35 a;* 
x"\ a7; + Sr — Carl te .. (2) 
a a)” 


such that the ®’s are universal functions of the boundary-layer 
coordinate, 7 = V (a1 v)y, for any given Prandtl Number. 
Since the Harvard Computational Laboratory has scheduled 
for the near future a relatively precise evaluation on the Mark I 
digital computer of some of these functions, the answers to the 
questions raised in the previous paragraph should be available 
soon. 

It seems appropriate to report, as a concluding remark, that 
we have been able more recently to include frictional heating 
effects in the “universal function”’ type of thermal boundary 
layer solution. 


’ 
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* Similar results for the so-called asymmetrical case, w=aix + aox? 


a3x* + ... have been obtained 


Interaction of Flow Pulsations and Flame 
Propagation{ 


G. H. Markstein 


Head, Combustion Section, Cornell Aeronautical Laboratory, Inc., 


Buffalo, N.Y. 
March 22, 1951 


I PRACTICALLY ALL TYPES of jet propulsion devices, pulsating 
flows occur under some operating conditions. The suppres 
sion of oscillations in steady-flow devices and the control of their 
amplitude in pulsating-flow engines have so far remained largely 
a matter of trial-and-error methods. The study of the phenome 
non in actual engines is complicated by the simultaneous occur 
rence of as yet poorly understood interactions between flow pul 
sations and diffuser operation, fuel injection, mixing, flame hold 
ing, and flame propagation. Although all contributing factors 
must be taken into account in a complete analysis, experimental 

t This research was conducted under the auspices of Project SQUID 


jointly sponsored by the Office of Naval Research and the Office of Air Re 
search under contract N6-ori-119, Task Order 1 


-JUNE, 19851 


and theoretical studies that attempt to separate the various as- 
pects appear to be indispensable before the final synthesis cay 
be carried out. 

The spontaneous excitation of organ-pipe oscillations by a 
flame propagating in a premixed combustible gas from the open 
to the closed end of a straight tube! seems to provide a cop. 
venient approach to the separate study of the interactions be. 
tween flow pulsations and flame propagation. The writer js 
currently studying vibratory flame motion in tubes with par. 
ticular attention to the hitherto neglected aspect of flame shape 
or structure. A detailed discussion of the apparatus, techniques, 
and results will be published elsewhere. For the present pur- 
pose, it suffices to mention the tollowing qualitative results, ob- 
tained by means of high-speed shadowgraph and schlieren cine- 
matography and by oscillographic recording of pressure tran- 
sients: 

(1) The vibratory phase of flame propagation was always 
accompanied by periodic appearance and disappearance of 
cellular flame structure.2~+ The positions of the centers of the 
cells appearing in any particular cycle coincided with the inter- 
sections of the cell boundaries in the preceding cycle. Thus, 
the motion of the flame front could be regarded as composed of 
standing waves with twice the period of the associated organ-pipe 
oscillation. The size of the cells increased with increasing ampli- 
tude of the vibration. 

(II) The transients of pressure measured at the closed end 
of a tube of 120-cm. length during flame propagation from the 
open toward the closed end were recorded for various combustible 
mixtures. These pressure records showed that for certain ranges 
of composition the oscillations started with a frequency corre- 
sponding to the first harmonic of the tube when the flame had 
traveled about one-third of the tube length from the top, while 
for the other compositions they started with the fundamental fre 
quency and only after the flame had passed through about haif 
the length. Methane-air mixtures showed no excitation of the 
first harmonic, while for higher hydrocarbons they appeared on 
the rich side of stoichiometric and for hydrogen-air and hydrogen 
hydrocarbon-air mixtures on the lean side. Similar trends were 
previously found for the appearance of cellular structure in slow 
burning mixtures.* 4 Detailed studies now under way indicate 
that the limits of composition for excitation of the first harmonic 
lie always within the limits of spontaneous cellular structure (i.e., 
cell structure present before the onset of oscillations ) 

These observations suggested a generalization of the first-order 
perturbation theory of flame front stability* to include the case 
of vibratory flame motion. The previous treatment yielded 


solutions of the form 
; ‘ 
exp (ify) exp (ol 


where the stability parameter 6(/) was determined by a char 
acteristic equation 
6? + 2B(h)6 + C(h) = 0 l 

Here, the coordinate y is parallel to the flame tront; 2 = 27/), 
where \ is the wave length of the disturbance in the direction y 
andi = VY —1. 

In view of the effect of acceleration on interface stability,® it was 
anticipated that the observed phenomena were caused by the 
alternating accelerations associated with the oscillation of the 
gas column. The introduction of a_ periodic acceleration 
—w°A cos wi acting normal to the flame front led to the follow- 


ing solution: 


£ = const. exp (hy): exp (—Bt)- V{(w/2)t,a,q] y 
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where 

F = deflection of the flame front normal to its equilib 
rium plane x = 0 

@ = angular frequency of the oscillation of the gas 
column 

A = amplitude of the oscillation of the gas column in 
the vicinity of the flame front 

B = defined by Eq. (1) 

P = ratio of densities of unburned and burned gases 

Y(z, a, g) = asolution of the Mathieu equation® 


(d?¥/dz?) + (a — 2q cos 22) Y = 0 


2 = (w/2)t 


(4) 


It is of particular interest in connection with the experimental 
results to note that, for certain regions of the (a, g) plane, this 
equation has unstable solutions of fundamental angular frequency 

9—i.e., of the form 


w/a , 


Y = const. exp (uz) . ( 


r 


+,exp [((2r+ 1)iz] (5 


with real positive p(a, g) : The condition for instability of solu 


tion (2) is then u(w/2) — B > 0. For any values of the ampli 
tude and frequency of the oscillation, Eqs. (8) determine the 
ranges of X = 22/h for which this condition is fulfilled Pre 
liminary computations have shown that for the observed ranges 
most cases < 1. 


of cell sizes and frequencies ja} is small—in 


Thus, its variation with \ is of minor significance. Taking as 
the most likely region of instability that between the character 
istic lines? b} and a; and a curve about halfway between these 
lines as that of maximum instability, one obtains (for small a) 
g~4. Withe 


magnitude with the observed relation between cell size and ampli 


= 5 this leads to A 2A, which agrees in order of 
tude. 

There remains to explain on the basis of the above analysis 
why the oscillations are spontaneously excited The condition 
for thermal excitation of organ-pipe oscillations has been given 
by Rayleigh:’ Heat 
along the tube with a component in phase with the pressure os 
The distortion £ of the flame tront of 


must be added periodically somewhere 
cillation at that point 
angular frequency w/2 is accompanied by changes of flame sur 
face area. These comprise both an increase of average area, 
responsible for the well-known increase of flame speed during 
vibratory flame motion, and a periodic change of area and, thus, 
of heat release with angular frequency w. The latter provides 
the feedback necessary for spontaneous oscillation; whether a 
certain frequency will be excited in any particular position of the 
flame front will depend on whether Rayleigh’s phase condition is 
fulfilled 

It is obvious that solutions for — of fundamental frequency w 
would give rise to heat-release changes of frequency 2 and thus 
could not lead to excited oscillations. This explains a posteriori 
why those regions of instability in the (a, g) plane corresponding 


to solutions 


1 


) = const exp (2r1 


eXp (us > ( 


must be The 
stable and unstable flames discussed under item (11) becomes at 


excluded. observed difference between initially 
least qualitatively understandable from the appearance of the 
stability parameter 6 in Eq. (3) for a. A detailed discussion of 


this analysis will be published at a later date. 
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Concerning a Definite Integral 


Yudell L. Luke and Dolores Ufford 
Midwest Research Institute, Kansas City, Mo. 
March 12, 1951 


I A PREVIOUS NOTE, Fettis! has shown that the definite integral 


priate e.= ile ll ie 
F(iz) = es \ . dl 


e™ sinh 6; 1 + sinh 6 — cosh @) cosh @ 
nie dé (1) 


0 sinh 6 


can be expressed in terms of tabulated functions except for one 
part that is expressed as a power series. It is the purpose of this 
note to show that Eq. (1) can be completely expressed in terms of 


tabulated functions. We also correct an error in his final re 


sult. For the sake of clarity, the notation of Fettis is employed 


The exceptional part referred to is 


“a | nr F\(x 


W\ Xx = dx 

J0 1 2 & 
x {* Hix . ; 

=> . dx = = Hy(x) — H,'(x)] dx 
*, x - Jo 

"% 

T 

= Hox) dx — Hy(x 2 
-“LJ0 


where //,,(x) is the Struve function It then follows that 


] 
F(x) = 1 — y — log 2x + + 
x 
T « ' 
‘ [Ho(x) — Yo(x)| dx + Yi(x) — Ai(x (3) 
21 Jo f 
where J,,(x) is the Bessel function of the second kind. An alter 
native expression for Eq. (3) obtains by noting that 
= - 
Vo(x) dx = xVo(x) + 5x Hix) ¥ilx) — Voix) Aix (4 
J 0 - 
Replacing x by 72, it is easy to show that 
Fr(1 = 1— y — log 2: + 
£ To - | dz +1 I(2 ) 
2t Jo f 
r 1 “ 
Fj(ts) = —- — + Ail + K | 6 


where J,,(x) and K,(x) are the modified Bessel functions of the 


first and second kinds, respectively, and where L,(x) is the modi 


fied Struve function. All functions appearing in Eqs. (3), (5 
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and (6) are tabulated.?~§ 
and (6) obtain by observing that 


z 
Is) dz = * s[folz)L—(s) — Loz)hi(2)] (7) 
0 - 
z ’ x : ; ; 
Koz) dz = 5 21Ko(2)L (sz) + Lo(s)Ki(z)] (8) 
0 « 


Fettis’ results for Fr(iz) and F;(iz) are incorrect, since he replaces 
E,\(iz) by L,(z). The correct relation is E,(iz) = (2/r) + Li(z) = 
L-:(2). 

Employing these results, the authors have tabulated values of 
Fr(iz) and F,(iz) and find that they agree with those given by 
Dingle and Kuessner.’ The latter tabulate S(z) = 1 — iz Fiz). 
The authors have also tabulated some auxiliary functions in con- 
nection with Eqs. (5) and (6) to facilitate interpolation in the 
vicinity of the origin. It is contemplated to deposit these data 
in the UMT file of the journal, Mathematical Tables and Other 
Aids to Computation. 
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On the Virtual Mass of Water Associated with 
an Immersing Wedge* 


John D, Pierson 

Hydrodynamics Development Engineer, The Glenn L. Martin Com- 
pany, Baltimore, Md.; Consultant to Experimental Towing Tank, 
Stevens Institute of Technology, Hoboken, N.J. 


March 28, 1951 


I A RECENT REPORT! on the penetration of a fluid surface by a 
wedge, the iterative solutions for the free surface shapes and 
impact forces were carried out for deadrise angles, 8, of 20°, 
30°, 40°, and 50°. These solutions were based on potential flow 
of an ideal fluid without gravity forces, but otherwise the actual 
flow conditions of spray formation and constant-pressure free 
surface were properly accounted for. During the course of the 
analysis, the potential distribution over the free surface and the 
solid boundary was calculated for use in the determination of 
pressures on the wedge. It now appears that an extension of the 
analysis to the determination of the fluid momentum would be 
of interest for comparison with the concepts of virtual water 
mass for float impacts. 

The penetration of the water surface by a solid body generates 
a fluid motion throughout the volume of water. Under the usual 

* The material reported herein was obtained in the course of the Research 


on Planing Surfaces at the Experimental Towing Tank of Stevens Institute 
of Technology under Contract with the Office of Naval Research 
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assumptions of the ideal fluid, initially at rest, this motion jg , 
potential flow described by the velocity potential, ¢. Th 
momentum of this field of flow is the summation over the entire 
region of the incremental masses of water, p d(vol.), times the 
velocity vector, V¢: 
Momentum of fluid = oS S Sve d (vol.) 
(vol.) 

It has been noted in a number of analytical developments for 
flows around flat plates, airfoils, etc., that a simple expression 
for the above integral may be found of the form 


mM.V = eS S Se d(vol.) (1 


(vol.) 

where V is the velocity of the body in the direction of the force 
exerted on the fluid and m, is the virtual mass of the fluid. [p 
other words, the integral of the momentum of individual par. 
ticles of the fluid flow is expressed in terms of a fictitious mass, 
My, Moving with the body. 

With the aid of Gauss’ theorem, the volume integral may be 
reduced to a more simple integral over the surface, 

SSS V¢dvol.) = —p JS Sed 
(vol.) S 
where df is a surface element as a vector (¢ df normal to surface), 
For the two-dimensional analysis, the surface integral is further 
reduced to the integral over a contour of unit width, 
m,».V = —p f oa (2 

The potential at any point in the fluid is a function of its rela- 
tive position, the size of the field (scale), and the velocity of pene- 
tration. Thus, 

@ = Vco(R, A) 

where c = the nominal wetted half width and (R, 6) = the co- 
ordinates of a position (relative to the width c). 

The surface element, df, is a function only of relative position 
Thus, 


df = c df (R, @) 


and the size of the field. 


The momentum can then be written in a form to separate the 
dimensional from the nondimensional parts, 


m.,V = —p f Vea(R, 0)c df (R, @) 
pve f 6’ df’ (3) 


where ¢’ and df’ are the functions of relative position only and 
the V, c terms may be removed from under the integral because 
they are scale functions and not dependent on the relative posi- 


tion. 


The virtual water mass per unit length is thus given by 
My = —pc* f o’ df’ (4) 


The values of the integral for various deadrise angles may be 
found by graphical integration of the potential distribution over 
the body and free surface given in reference 1. These calculated 
values are shown in Fig. 1 where similar values of m,,./(a/2)pc? 
from earlier methods are shown for comparison 


von Karman? Mwy = (x/2)p (2c/xr)? 
Wagner? Mw = (w/2)pc? 
Kreps? My = (w/2)pc?{1 — (B/r)| 


: : *e Jf 2/2 tan B\? 
Wagner-Sydow?) my = = pc? — | 
7 2 2 7 
z e Jf 2/2 tan B\? 
Mayo? Mm.» = 0.82 — pc? — 1 
2 2 T 


Monaghan? My = (w/2)pc?|1 — (B/r) |} 
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Fic. 1. Variation of virtual water mass with deadrise angle 
It will be noted that the calculated values of the virtual mass 
(Table 1), based upon the surface potential distribution, lie be- 
tween the approximate formulas proposed by Mayo, Wagner, 
and Monaghan. At low deadrise angles the trend is toward the 
former, while in the high deadrise-angle range the calculated 
values lie closer to the latter two. Monaghan’s analysis’ is based 
upon an analytical study of the flow about a submerged diamond- 
shaped object (which neglects the formation of spray) and, thus, 
will apply only to the extent that such a flow is representative of 
the flow for the immersing wedge. The approximate variation 
of impact force with deadrise angle proposed by Wagner? (inter- 
preted in terms of virtual mass by Sydow) is based upon a de- 
tailed calculation of the flow for one angle of deadrise (18°) and 
a fairing into the conditions at 0° and 90° deadrise. (In the 
present analysis, the virtual mass has been computed directly 
from the flow field for four deadrise angles.) This Wagner- 
Sydow determination of virtual mass is reduced by a constant 
empirical factor (0.82) determined by Mayo from correlation with 
planing data. Later study’ has indicated that the magnitude 
of this factor to be applied to the Wagner-Sydow relation for the 
virtual mass may be much closer to 1 at higher deadrise angles 
Thus it appears that the calculated values of the virtual mass 
based upon the detailed analysis of the flow field (as represented 
by the surface potential distribution) are not only consistent 
with the theoretical flow conditions for the immersing wedge but 
also agree closely‘with the approximate empirical values that have 


been obtained 


TABLE 1 





von Wagner- Present 
8 Karman Wagner Kreps Sydow Mayo Monaghan Analysis 
10 0.405 1.00 0.945 0.802 0.657 0.844 
20 0.405 1.00 0.89 0.66 0.54 0.703 0.575 
30 0.405 1.00 0.823 0.541 0.444 0.58 0.483 
40 0.405 1.00 0.778 0.446 0.366 0.471 0.40 
0 0.405 1.00 0.723 0.367 0.300 0.376 0.34 
60 0.405 1.00 0.667 0.304 0.25 0.297 
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The Source-Flow in Supersonic and 
Incompressible Theories 


H. Serbin 
Professor, School of Aeronautics, Purdue University, Lafayette, Ind. 
March 26, 1951 
I THE THEORY OF A PERFECT FLUID in three dimensions, one 

defines a source of strength Q at a point P(0,0,0) as: (a) a 
potential flow depending only on the distance 

y = (x? + y*? + 2%) 

such that the (volumetric) rate of flow across a sphere with center 
at P is equal toQ. 

From this, it follows that (b) a source of strength Q is defined 
by the potential 


¢ = —Q/(4rr) (1) 


From definition (a), one infers that a continuous distribution of 
sources, on a plane, of intensity g will produce a flow such that the 
normal velocity at the plane is g/2 outward. Therefore, (A) if 
the plane is immersed in a strong uniform flow field of velocity 
U(U > q) so as to be parallel to the stream lines, the stream 
lines are deflected by the presence of the sources so that the re- 
sulting slope is g/(2U). 

In a similar way, consideration of a continuous distribution of 
weak sources g along a straight line immersed in a uniform flow 
so as to be parallel to the stream lines shows that (B) the flow 
produced by a continuous distribution of sources of intensity 
q(q < U) is such as to alter the stream lines at small radius r 
to the slope g/(2r LU’). 

Statements (A) and (B) are important for the thin airfoil 
theory and body theory, respectively. It is worth while to con- 
sider the extension of these results to linearized supersonic flow. 

Write the potential equation in the form 


Grr — @ ? (Gyy + ox) = 0 (a? = M*? — 1) (2) 


and define s by 


[x? — a%(y? + 2?)]’ 


5= 


It does not appear feasible to start from a definition of a source 
comparable to (a); instead, it is a common practice to define a 
source, analogous to (b), by the solution of Eq. (2) 

¢ = A/s (3) 
In this form, the numerical coefficient —Q/(47) has been replaced 
by a new constant A, but there is no physical significance at 
tached to this constant. 


On the other hand, suppose one defines a source Q in the stream 


by 
¢g = —Q/(2ns) (4) 
where we use the number 2 instead of 4 as in Eq. (1). Then 


statement (A) holds, as has been shown by Puckett.' This sug- 


gests a parallelism between the results of the perfect flow and the 
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linearized supersonic theory (this parallelism appears to be the 
basis of Hadamard’s work) 

However, an examination of the von Karman-Moore theory 
of rotationally symmetric linearized supersonic flow reveals a 
difficulty with sign. 
stream lines at radius r is given by 


According to their results, the slope of the 


dr/dx = —q(x)/(2rrU) (5) 
which differs in sign from result (B) and where r = (y? + g?)'/2, 


The derivation leading to Eq. (5) is based upon the following 
sequence of steps :? 


t—~ar 
g(x,7r) = f f(&) (xe — &)? — ar?]-'/ dé 
0 


(f(t) = — g(é)/2r] (6) 


[(x — —&)/(ar) =coshu] (7) 


= S fx — ar cosh u) du, 


=—-—a St f’'(x — ar cosh 1) cosh u du (8) 


1 x—ar 
[ f'(E) (x — &) [(x% — £)? — a?r?] "/2 dé (9) 
rJo 
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where the integration in u is from the lower limit cosh™! (x/ar) 
to the upper limit 0. The transformation £ — w is two-valued, 
however, and there is an ambiguity attached to the sequence of 


Eqs. (7) and (8). If one chooses the branch uw S 0, then 


( +\9 9.91!/o . 
[(x — £)? — a®r?]/? = —ar sinh u 


d& = —ar sinh u du 


du = [(x E)2 — qr?]—'/2 dé 


The choice of the negative branch is consistent with Eqs. (6)-(8) 
but not with Eq. (9), which, according to our analysis, should 
have a negative sign. With this modification, statement (B) 
holds for linearized supersonic flow, 

The change in sign does not affect the wave drag, which is of 
quadratic character. However, because the discrepancy in sign 
mars the theory, it was felt advisable to draw it to the attention 
of your readers, 
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Symposium on Fluid Dynamics 


The Fourth Symposium on Applied Mathematics of the American Mathematical Society 
(cosponsored by the University of Maryland and the U.S. Naval Ordnance Laboratory) will be 
held at the University of Maryland, College Park, Md., and the U.S.N.O.L., White Oak, Md., 


The general topic of the Symposium is “Fluid Dynamics,’’ and sessions are planned on 
Turbulence, Foundations, Compressible and Incompressible Flow. 

I.A.S. members interested in attending may secure additional information by writing to 
M. H. Martin, University of Maryland, College Park, Md. 














